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ABSTRACT 


This  research  modifies  the  existing  finite  element  formulation  of 
a  potential  energy  based  large  deformation  and  moderate  rotation  theory. 
Herroitian  shape  functions  replace  the  existing  linear  bending  angle 
interpolations.  Negligible  differences  between  the  two  formulations 
indicate  the  underlying  kinematics  limit  the  accuracy,  not  the  finite 
element  interpolations.  Using  the  new  program,  numerous  nonlinear  arch 
geometries  are  modeled  to  investigate  the  effects  of  arc  length  and 
thickness  variations .  Local  and  global  snapping  phenomena  are  captured 
as  well  as  through  the  thickness  shear  driven  nonlinearities. 
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INVESTIGATION  INTO  THE  BEHAVIOR  OF  GEOMETRICALLY 


NONLINEAR  COMPOSITE  ARCHES 


I .  Introduction 


1 ■ 1  Background 

Composite  shells,  plates,  arches,  and  beams  are  the  building 
blocks  for  the  majority  of  today's  complex  aerospace  structures.  The 
critical  nature  of  such  structures  requires  accurate  prediction  of  their 
behavior  under  any  loading  conditions.  However,  analysis  beyond  the 
basic  linear  geometric  and  material  regions  is  complex  and  time 
consiuning . 

One  type  of  complex  shell  analysis  combines  small  strains  with 
large  displacements  and  rotations.  In  this  situation,  the  material 
remains  linearly  elastic  while  the  geometry  becomes  non-linear. 

Geometric  non-linearity  causes  dramatic  structural  stiffness  changes 
without  significant  material  changes.  Palazotto  and  Dennis  present  the 
Simplified  Large  Displacement /Rotation  (SLR)  theory  and  complementary 
FORTRAN  code  to  analyze  shells  experiencing  these  loading  conditions 
(19) .  SLR  theory  determines  the  equilibrium  path  of  orthotropic  shells 
using  a  total  Lagrangian  approach  and  includes  a  parabolic  transverse 
shear  stress  distribution.  Equations  are  solved  using  the  finite 
element  method.  The  result  is  a  FORTRAN  code  capable  of  analyzing 
linearly  elastic  behavior  of  isotropic  and  laminated  orthotropic 
cylindrical  shell  structures  undergoing  geometric  non-linearity  with 
large  deformations  and  moderate  rotations. 

Creaghan  reduced  SLR  theory  to  one  dimension  for  the  analysis  of 
beams  and  arches  (3;  4).  The  resulting  FORTRAN  code  combines  all  the 
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features  of  Palazotto  and  Dennis's  work  with  extensible  mid-plane  strain 
terms  from  Smith  (25;  26) .  Creaghan  also  include  displacement  control 
and  Riks  methods  for  solving  non-linear  equations.  These  modifications 
allow  uninhibited  observation  of  SLR  displacement  and  rotation  limits. 
Creaghan 's  theory  produces  accurate  solutions  through  23  degrees  of 
rotation. 

Miller  improved  Creaghan 's  rotation  limits  to  approximately  45 
degrees  (13).  Miller  replaced  Creaghan’s  small  angle  approximation  for 
the  bending  angle  (\|/)  with  a  truncated  series  representation  of  the 
tangent  function,  thereby  introducing  large  rotation  kinematic  theory  to 
Creaghan ' s  formulation . 

The  purpose  of  this  thesis  is  to  improve  the  rotation  limit  of 
Miller's  code.  Miller's  adaptation  of  the  finite  element  method  (FEM) 
will  be  altered  to  include  higher  order  shape  functions  for  \(f. 

Refining  the  finite  element  model  will  improve  solutions  within  the 
limits  of  Miller's  theory.  The  resulting  FORTRAN  code  will  then  be  used 
to  analyze  arch  classification  parameters. 

1 . 2  Preview 

Thickness  and  depth  primarily  determine  the  behavior  of  shells  and 
arches.  For  example,  shear  plays  a  large  role  in  thick  structures  but 
is  negligible  in  thin  shells;  therefore,  a  theory  which  considers 
transverse  shear  should  be  used  for  analysis.  Also,  it  is  important  to 
note  a  shallow  shell /arch  behaves  very  much  like  a  plate/beam  whereas  a 
deep  shell/arch  does  not.  Obviously,  the  geometry  of  the  structure 
should  be  classified  before  analysis  begins. 

Huang  proposes  thin  arches  have  a  very  small  physical  thickness 
{h)  vs.  radius  of  curvature  (r)  ratio  (h/r  «1)  (10).  Miller  also 
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considers  the  half  chord  length  (c) ,  but  his  decision  value  {c/h<25)  is 
arbitrary  (13). 

Smith  proposes  a  shallow  shell  has  a  small  depth  (8)  vs.  half 
chord  length  ratio  (5/c<0.3)  (25).  Fung  and  Kaplan  propose  X.  =  2[3(l-u2)]-25('5 

/h)-^  as  a  classification  criteria,  limited  to  small  opening  angle  (a«l) 
(8) .  Kaplan  revises  this  definition  to  remove  the  opening  angle 
limitation  resulting  in  X  =  (12).  Miller  uses  A,  >8  to 

define  deep  shells  (13) . 

The  complex  nature  of  large  deformation  and  rotation  theory 
requires  solving  niimerous  simultaneous,  non-linear,  differential 
equations.  Even  with  simplifying  assumptions,  these  equations  remain 
complex  and  extensive.  For  this  reason,  most  theories  employ  FEM  for 
solution. 

Accurate  FEM  solutions  require  proper  shape  functions.  Cook, 
Gallagher,  and  Zienkiewicz  all  discuss  shape  functions  at  great  length 
(2;  9;  29).  Shape  functions  interpolate  quantities  within  an  element 
and  are  generally  classified  as  C®  or  continuous.  All  shape 
functions  must  be  continuous;  that  is,  the  quantity  they  interpolate 
must  be  continuous  at  each  node.  If  a  quantity  is  represented  by 
shape  functions,  the  first  derivative  of  the  quantity  is  continuous  as 
well.  Obviously,  or  Hermitian  shape  functions  are  required  if  a 
quantity  and  its  first  derivative  must  have  interelement  continuity. 
Miller  uses  linear  shape  functions  to  estimate  \\i  (13)  .  As  such,  the 
first  derivative  of  \|/  (\)^2)  is  constant  within  a  given  element  and 
discontinuous  at  each  node.  The  current  work  will  employ  both  quadratic 
and  Hermitian  interpolations  of  \|/. 
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1.3  Literature  Review 


Although  primarily  concerned  with  non-linear  shell  theory,  this 
thesis  also  addresses  arch  classification  parameters  and  finite  element 
method  shape  functions.  Therefore,  a  brief  review  of  these  topics  is  in 
order . 

Kirchhoff's  theory  is  one  of  the  first  for  analysis  of  thin  plates 
(22) .  His  assumptions  allow  two  dimensional  analysis  of  three 
dimensional  plate  behavior.  He  assumed  an  inextensible  mid  plane  and 
negligible  normal  strains.  Kirchhoff  also  eliminated  transverse  shear 
effects  by  disallowing  cross  sectional  warping. 

Love  took  the  next  logical  step  in  the  development  of  shell  theory 
(22).  He  adapted  Kirchhoff's  plate  theory  to  curved  shells.  The 
Kirchhoff -Love  theory  obviously  neglects  transverse  shear  strain 
effects.  As  a  result,  both  Kirchhoff's  theory  and  Love's  revision  only 
apply  to  thin,  isotropic  plates  experiencing  small  deflections. 

Reissner  and  Mindlin  improved  Kirchhoff -Love  theory  through  the 
addition  of  transverse  shear  considerations  (14;  21) .  Reissner-Mindlin 
(RM)  theory  assumes  constant  through  the  thickness  shear  and,  as  a 
result,  cross  sections  normal  to  the  datum  surface  rotate  but  do  not 
warp.  This  first  order  theory  can  cause  shear  locking  in  a  finite 
element  formulation  because  boundary  conditions  are  not  satisfied  on  the 
top  and  bottom  of  the  shell.  Shear  correction  factors  minimize  this 
condition. 

Reddy  adds  a  parabolic  transverse  shear  strain  representation  to 
RM  theory  to  eliminate  shear  locking  (20) .  Palazotto  and  Dennis, 
Creaghan,  and  Miller  employ  Reddy's  theory  (3;  13;  19).  As  such,  it  is 
retained  here . 
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Beam  and  arch  theory  spawn  an  enormous  amount  of  research.  As 
with  shell  theory,  arch  theory  refinements  are  directed  towards 
improving  displacement  and  rotation  limits.  As  a  baseline  for 
comparison,  this  thesis  retains  Miller's  theory  without  changes  (13). 

His  theory  combines  total  Lagrangian  kinematics  with  the  aforementioned 
trigonometric  representation  of  the  bending  angle  and  vector 
decomposition  of  the  displacement  equations.  Miller  also  includes 
parabolic  transverse  shear,  mid-plane  extensibility,  and  all  non-linear 
in-plane  strain  terms.  A  finite  element  formulation  approximates 
solutions . 

Oliver  and  Onate  present  a  total  Lagrangian  finite  element 
formulation  which  accurately  calculates  large  displacements  and 
rotations  (18) .  They  include  a  parabolic  shear  representation  and  mid 
plane  extensibility.  Degenerated  shell  elements  similar  to  SLR  theory 
are  used  as  well  as  exact  trigonometric  relations .  This  is  a  stress 
based  theory  whereas  SLR  is  strain  based. 

Surana  also  uses  a  vector  decomposition  of  the  displacement 
equations  (27)  .  Total  Lagrangian  kinematics  and  a  finite  element 
representation  are  employed  as  well . 

Noguchi  and  Hisada  use  a  finite  element  model  to  analyze  post 
buckling  shell  response  (17)  .  Their  vector  based,  total  Lagrangian 
kinematics  provide  a  good  comparison  to  the  current  derivation. 

Mondkar  and  Powell  develop  equations  of  motion  based  on  total 
Lagrangian  kinematics  (16)  .  They  retain  higher  order  terms  and  consider 
mid-plane  extension  and  transverse  shear.  Solutions  for  isotropic 
structures  are  found  using  FEM. 

Huddleston  uses  a  total  Lagrangian  formulation  with  no  small  angle 
ass\amptions  to  analyze  deep,  isotropic,  circular  arches  (11) . 
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Transverse  shear  is  neglected  while  mid  plane  extension  is  allowed. 
Huddleston  solves  simultaneous,  non-linear,  first-order  differential 
equations  rather  than  using  finite  element  or  finite  difference  methods. 

Sabir  and  Lock  also  neglect  transverse  shear  and  eliminate  some 
higher  order  mid-plane  strain  terms  (23).  Total  Lagrangian  kinematics 
are  used  in  conjunction  with  a  finite  element  solution  method  to  analyze 
isotropic,  circular  arches. 

DeDeppo  and  Schmidt  present  numerous  approaches  to  solving  non¬ 
linear  circular  arch  problems  (5;  6;  24),.  When  they  consider 
transverse  shear,  a  first  order  representation  is  used.  When  mid  plane 
stretching  is  considered,  higher  order  terms  are  neglected.  They  use 
total  Lagrangian  kinematics  with  no  small  angle  approximations  and 
finite  difference  solution  method. 

Epstein  and  Murray  retain  quadratic  in-plane  strain  terms  while 
neglecting  transverse  shear  strain  (7) .  Mid-plane  extension  is 
considered  with  higher  order  terms  retained.  Epstein  and  Murray  use 
finite  elements  and  a  total  Lagrangian  approach  to  examine  only 
isotropic  beams . 

Belytschko  and  Glaum  use  small  angle  approximations  in  an  updated 
Lagrangian  basis  (1) .  Mid-plane  stretch  is  allowed  but  higher  order 
terms  are  eliminated.  Also,  they  ignore  transverse  shear. 

Minguet  and  Dugundji  also  use  the  co-rotational  method;  however, 
Euler  angles  are  used  for  exact  rigid  body  kinematics  (15) .  Transverse 
shear  and  mid-plane  extensibility  are  considered  without  higher  order 
terms . 
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II .  Theory 


2 ■ 1  Overview 

The  current  work  is  based  on  Miller's  kinematic  development, 
constitutive  relations,  strain  relations,  and  beam  energy  considerations 
(13) .  His  work  is  reviewed  here  for  completeness  and  to  highlight 
certain  inaccuracies  in  the  theory.  Miller's  FEM  formulation  is 
modified  to  include  a  higher  order  shape  function  for  the  bending  angle. 
Numerical  solution  methods  and  arch  geometry  definitions  are  discussed 
as  well. 

2 . 2  Kinematics 

Miller  derives  arch  displacement  equations  in  a  vector  format. 

The  first  step  considers  displacements  due  to  pure  bending  as  shown  in 
Figure  2-1.  Points  j  and  k  are  on  the  outer  surface  of  the  beam  and 


Figure  2-1  Flat  Beam  Bending  Rotation 


point  i  is  on  the  mid  plane.  The  V  coordinates  are  initially  parallel 
to  the  original  system  and  move  with  the  normal  of  the  cross  section 
during  deformation.  The  e  coordinates  are  the  global  coordinates  of  the 
total  Lagrangian  system.  The  vector  remains  parallel  to  the  cross 


2-1 


sectional  thickness  during  bending.  The  magnitude  of  Vji  equals  the  beam 
thickness  {h)  .  Yy,  remains  perpendicular  to  Vai  during  pure  bending  and  a 
is  the  rotation  angle  about  Vn.  The  natural  coordinate  ^  is  aligned  in 
the  63  direction  with  magnitude  1  at  the  top  of  the  beam  (positive  63 
direction)  and  -1  at  the  bottom.  is  expressed  in  global  coordinates 

as : 

Ki=(3',-3’*K-(Zy-Z*)^3  <2-l) 


with  direction  cosines: 


_  1 

'yj-y^' 

."3/. 

h 

Eqn(2-2)  is  substituted  into  Eqn{2-1),  yielding: 


» 


(2-2) 


Vy"  =h(m"e2+nyej) 


(2-3) 


The  n  superscript  denotes  any  deformed  state.  The  position  of  any  point 
P  on  the  cross  section  is  described  by  its  location  from  the  beam  center 
line  and  the  location  of  point  i  on  the  mid  plane  as 


(2-4) 


The  total  Lagrangian  displacement  of  any  point  P  is: 


u(z)  = 


I 

1 _ 

1 

e 

- 1 

_ 1 

L^.-J 

(2-5) 


The  superscript  0  denotes  the  original,  undeformed  state.  Miller 
assumes  vertical  displacement  (w)  is  constant  through  the  thickness. 

As  shown  in  Figure  2-2,  P  only  displaces  horizontally  due  to  bending. 
This  assumption  implies  stretching  of  the  beam  normal.  Resulting  no2rmal 
strains  are  neglected.  Horizontal  displacement  is  assumed  in  order  to 
keep  w  constant  through  the  thickness .  The  resulting  strains  become 
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Figure  2-2  Bending  Motion 


significant  as  a  approaches  7t/2  .  This  assumption  changes  Eqn(2-3)  as 
follows : 


_  h  K,' 

cos(a,)L«3/. 

where 

ffiji  =-sin(af) 

=cos(a,.) 

Egn(2-6)  is  substituted  into  Eqn(2-5)  at  a(|=0  yielding: 


so  that 


if  r  M  r  T”  r  n® 

u  =  u.+—n -  - 

2  [cos(ai)[%J  Ln3,_ 

1  -tanfa,)!  „  „ 


(2-6) 


(2-7) 


where  m,.  represents  the  motion  of  point  P  due  to  the  motion  of  point  i 


on  the  midplane  and  [R]  is  the  bending  rotation  tensor.  The  general 


displacement  equation  using  natural  coordinates  is 


where  [1]  is  the  2X2  identity  matrix. 


(2-8) 
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Adapting  this  derivation  to  curved  beams  requires  consideration  of 
the  additional  motion  due  to  curvature.  Orthogonal  curvilinear 
coordinates  are  used  for  arches  as  shown  in  figure  2-3.  The  s  direction 


Figure  2-3  Global  Arch  Coordinates 

runs  along  the  arch  midplane  and  the  3  direction  points  toward  the 
center  of  curvature. 
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The  u,  term  in  Eqn  2-8  represents  the  displacement  of  an  arbitrary- 
point  caused  by  non-bending  motion.  For  straight  beams,  mid-plane 
motion  directly  correlates  to  Uj  for  any  arbitrary  point.  As  shown  in 
figure  2-4,  the  becomes  thickness  dependent  in  a  curvlinear  system 
therefore  requiring  metric  tensor  coefficients  (19;  22) .  In  this  case. 


the  motion  of  point  P  due  to  midplane  motion  is; 


M,  == 


v(l--) 

r 

w 


(2-9) 


where  v  and  w  are  the  motion  of  point  i  on  the  midplane.  Using  the 
shell  scale  factors  as  did  Palazotto  and  Dennis  yields  similar  results 
(19)  . 

Figure  2-5  shows  the  relationship  between  \|/  and  slope  (  w.j )  .  A 
positive  \)/  is  defined  as  the  angle  which  causes  a  point  with  a  positive 
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z 

Figure  2-5  Pure  Bending  Convention 


z  coordinate  to  move  in  the  positive  s  direction.  Since  a  is  defined 
using  the  right  hand  rule  in  a  positive  direction  about  the  1  axis, 


\|/  =  -a. 

Expressing  Eqn(2-8)  and  Eqn{2-9)  in  global  coordinates  yields 
these  displacement  components : 

«,  =  vf  1-—  |+ztan(\|/) 

'  I  rj  ^  (2-10) 

Mj  =  W 

These  are  the  displacement  relationships  for  bending  only.  Instead  of 
using  a  small  angle  approximation  for  tan(\|/).  Miller  used  shell 
displacement  relationships  developed  by  Palazotto  and  Dennis  to  retain 
the  tangent  function  in  the  Ui  component (13 ;  19).  As  a  result,  the  in 
plane  displacement  of  Eqn(2-10)  becomes: 

M2(j,z)  =  v|^l-— j+ztan(\|/)  +  z’A:(\j/  +w,2 ) 

where:  (2-11) 
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similar  to  \|/,  a  positive  shear  angle  P  causes  a  point  with  a 
positive  z  coordinate  to  move  in  the  positive  s  direction.  Figure  2-6 
shows  the  relationship  between  P  and  w,2.  Mathematically,  w,2,  P,  and 
are  related  by: 

w,j+v=-p  (2-12) 

Eqn{2-12)  is  used  to  find  P  after  and  \);  are  found  using  FEM.  The 
sign  convention  of  these  three  quantities  is  shown  in  figure  2-7. 


Figure  2-6  Angular  Shear  convention 


The  tangent  function  in  Eqn{2-ll)  is  approximated  using  the 
following  series  expansion: 

tan(v )  =  X -l'~'  ~ (2-13) 

(=0  "  • 

where  the  Bn's  are  Bernoulli  numbers  and  \|/  is  in  radians  (12)  .  This 
series  expansion  is  truncated  at  two  terms  yielding: 


tan(\|/)  =  \i/+-\|/^ 


(2-14) 


2-7 


Figure  2-7  Slope,  Shear,  and  Bending  Relationship 


Substituting  Eqn(2-14)  into  Eqn(2-ll)  for  M2,  one  obtains  the  final 
displacement  components : 

M[  =  0 

v^l— ^j  +  z(\|/  +-i\j/^)  +  Z^A:(\|/  +w,2 )  (2-15) 

Mj (j) =  W 

where  v  and  w  are  midplane  displacements  in  the  2  and  3  directions, 
respectively,  z  is  the  distance  from  the  midplane, \j/  is  angle  due  to 
bending,  and  w,2  is  the  slope  of  the  midplane  tangent. 
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2.3  Constitutive  Relations 


Table  2.1  presents  the  explicit  and  contracted  notation  used  for 
both  the  stress  and  strain  tensors. 

A  fiber  composite  lamina  has  the  local  coordinate  system  defined 
as  shown  in  figure  2-8.  The  1'  direction  is  aligned  with  the  fibers  of 


Table  2-1  Contraction  Definitions 


Stress 

Strain 

Contracted 

Explicit 

Contracted 

Explicit 

a, 

0,1 

e. 

£11 

02 

CT22 

El 

£22 

03 

033 

£3 

£33 

04 

023 

£4 

2623 

05 

0,3 

Es 

2£i3 

06 

0,2 

Et 

2£i2 

the  composite.  The  2'  and  3'  directions  are  the  transverse  direction. 

In  this  work,  the  3  and  3'  directions  are  aligned  and  point  into  the 
page. 

Following  these  conventions,  the  stress  and  strain  for  a  composite 
lamina  are; 


'Qn 

Qi2 

0 

0 

0  ■ 

<^2 

Q\2 

Q22 

0 

0 

0 

£2 

^6 

.= 

0 

0 

Qee 

0 

0 

£5 

<^4 

0 

0 

0 

644 

0 

£4 

0 

0 

0 

0 

to 

1 _ 

(2-161 


The  Qij's  represent  the  ply  stiffness  expressed  in  material  coordinates, 
the  a' ' s  are  ply  stresses,  and  the  ' s  are  the  strains. 

These  relations  assume  plane  stress  conditions  exist,  therefore  is 
negligible.  Through  the  thickness  strain  (63)  is  found  through 
constitutive  relations  as  shown  by  Palazotto  and  Dennis  (19).  The  Qij's 
in  terms  of  engineering  constants  are; 

E, _ 

V..  1-V..V,,  l-v„v,. 

(2-17) 


Q  _  ~1  Q  ^  ^12^2  -  ^21^1 

1— VijVji  1~'^I2'’21 


G22  ~ 


l-V,2Vj, 


Qu  ~  ^23  Q55  ~  ^13  Qee  ~  ^12 


In  these  relations,  E^‘ s  are  Young's  moduli,  Gi j ’ s  are  shear  moduli, 
and  V;;'s  are  Poisson's  ratios.  All  values  are  expressed  in  material 

V 

coordinates . 

When  constructing  a  laminate  composite,  the  stiffness  (2y)  of 
each  lamina  must  be  converted  to  global  coordinates.  This 
transformation  is  accomplished  using  the  angle  6  defined  in  figure  2-8 
and  the  transformation  equations  found  in  Reddy  (20) .  These  equations 
are  presented  in  Appendix  A  and  give  the  following  constitutive 
relations  for  the  ^th  ply; 
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(2-18) 


£2 

■£6  ■ 

£4 

where  Qy  ' s  represent  the  plane  stress  stiffnesses  in  global 
coordinates . 

The  next  developmental  step  requires  simplifying  the  2-D  Eqn(2-18) 
to  1-D.  Because  this  work  examines  beams  and  arches  which  are  narrow  in 
the  1  direction,  normal  stresses  in  this  direction  ffi  are  negligible. 
Also,  because  the  beam/arch  is  narrow  and  traction  free  at  the  sides, 
the  in  plane  shear  stress  Oe  is  neglected.  Stress  and  strain  in  the  5 
direction  (Os  and  £5)  are  neglected  because  beam  twisting  is  not 
considered.  These  assumptions  reduce  Eqn(2-18)  to  the  following  1-D 
constitutive  relations  (3) : 


'a, 

Qn 

6i6 

0 

0 

^2 

Qi2 

Q22 

G26 

0 

0 

<^6 

-  = 

Qi6 

^26 

^66 

0 

0 

<^4 

0 

0 

0 

Qu 

0 

.<^5. 

k 

0 

0 

0 

0 

Q55. 

o 

<J2*=4*£2t 


where 


(2-19) 


Qlk  ~ 


Qil  ^ 

V/22  py 
V  v^ii , 


Because  twisted  and  torsional  stiffness  are  not  considered,  Eqn(2- 
19)  is  limited  to  balanced,  symmetric  lay-ups.  This  limitation  drives 

j2i6  to  zero  for  the  entire  laminate. 
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2.4  Strain  Relations 


Miller  and  Creaghan  both  employ  the  following  relationship  between 
Green  strain  and  shell  physical  strains  developed  by  Palazotto  and 
Dennis  (19) . 


e 


ij 


h. 

h,hj 


(2-20) 


In  Eqn(2-23),  the  Yy's  are  Green  strain  tensor  components  and  the  /i,'s 
are  shell  scale  factors.  The  previously  discussed  thiclcness  dependence 
of  Uj  drives  the  requirement  for  shell  scale  factors .  Because  only  in 
plane  normal  and  through  the  thiclcness  shear  strains  are  considered, 
only  the  two  Green  strain  tensor  terms  related  to  these  strains  are 
required. 


,  /ijM,  ,  , 

Y  22  —  ^2^2,2  ,  ^2,3  ,  ^2.1 


h, 


1 

+  — 

2 


1 

+— 

2 


^  «3  .  «1  r  ^ 

“2,2  ,  ^2,3  .  ^2,1 

hj  h,  j 


U. 


“3,2  ,  ^2,3 

V  *3 


V  1 
+  — 
2 


^  u 

^2  j 

^1,2  ~ ,  ^2,1 
hi 


J 


(2-21) 


Y  23  ~  2  ^^3^3.2  "^^2**2,3  **2^2,3  **3^3,2) 


1 

+  — 
2 


1 

+— 

2 


1 

+— 

2 


^  “3 ,  Y  “37  “1 7 

“2,3  ,  ^2,3  **2.2  7  ^2,3  7  ^2,1 

K  jl  h, 


«7 


**3,2  7  ^2,3 

V  ^3 


**2  7  **1  7 

**3,3  7  ^3,2  7  ^3,1 


/*. 


..  V 
2 


**1,2  7  ^2,1 
\  hi  J 


—!-Lu 

**1,3  7  "3,1 


(2-22) 


Here,  the  w,.  's  come  from  Eqn(2-15)  . 


2-12 


The  scale  factors  for  a  cylinder  with  curvature  in  only  the  2  or  s 


direction  are: 


h^=h,=l 

/ij  =1-- 


(2-23) 


Since  beam  width  does  not  affect  these  scale  factors,  they  also  apply  to 
a  curved  beam  or  arch.  These  factors  match  Eqn{2-9) ,  which  were 
developed  graphically  from  figure  2-4. 

Inserting  Eqn(2-15) ,  Eqn(2-21),  and  Eqn(2-23)  into  Eqn(2-20) 
yields  the  in  plane  normal  strain  £2  terms.  All  terms  of  Eqn(2-21)  are 
retained  because  in-plane  strain  dominates  thin  beam  bending  problems. 
Representing  the  scale  factors  with  truncated  Taylor  series  simplifies 
expansion  of  Eqn(2-20) .  The  resulting  sixty  expressions  are  shown  in 
Palazotto  and  Dennis  (19).  For  this  worlc,  only  the  following  nine  of 
the  sixty  scale  factor  approximations  are  required: 


1  , 

—  “  1  +  zc 
h. 


‘2,3  2 

■  a—c  —  C  Z 


"2,3 

hlh. 


~  -c-2c^z 


"2,3 

/lj/13 


•  c^  +c^z 


±.U2. 


•  +2c^z 


(2-24) 


'*2.3  2  .  3 

- ~c  +c  z 


h 


2,3  r,  2 

—  =  -c-2c  z 


where  c  =  1/r.  In  plane  normal  strain  is  expressed  in  terms  of  the  mid 
plane  strain  (8°)  and  functions  independent  of  Z  {X2i^  ' 


£2 


2p 


(2-25) 


The  mid  plane  strain  and  Xji'®  defined  as: 


e®  =  V2  -  wc+.5{v 2^  +'^,1  +v^c^  +  w^ c^)+vw 2C-V 2^0 
X2I  =V,2  +W2^C  +  W^C^  -C^(V2W-Wi'2)  +  V\|/C^  +^,2^,2 

-  c(v  _2  W  -\t/w_2 )  +V|/  ,2V '  +  v,2¥  .2¥ "  -  w-cii/  .2V " 

X22  =V,2C+-5('K,2^  +\(f  ^c^)  +  v\|/c^  -2c^(\|/  2M'-Vvv2)+V|/,2V2C-1.5(c‘'v^  h-c^Vj^) 

-  c  ^  V  2  +  2c  ^  (v  2  w  -  vw  2 )  +  c(v|/  2V  ^  +¥  .2^  ^  V  2 )  +  V]/  ,2 V  ^ 

X23  =*(^,22  +¥,2)  +  C\1/,2^  +\t/^C^  +^.2^22  +V  .2  )  +  (w^  +11^  )  “  H'fc(W,22  +^,2) 

+  W2^C(W2  +\|/)  +  c^v^  +c^V2^  -2c^(c^v\i;+V2M>_2) 

X24  =kc(W22  +\\f^2)  +  vkc\w2  +\^)+2kc\-WW22  -W,2V  +^,2^  +  ,2  ) 

+  *\|/  2(^22  +V,2)+V^C^(W2  +\ir)  +  V2fc(M',22  +V,2)  +  ^.2¥^(1  +  W,22) 

X25  =2fe[\/  2(W,22  +V,2)+VC^(W2  +V)-‘^.2(w',22  +V  ,2  )  “  2  +¥)+¥  ,2V  ^^,22  1 

Z26  =yK2'  +2vt',22¥.2  a  +c\w2^  +lw_2\f 
X27  =A:'c[(w22+V,2)'+c'‘(w2  +V|f)"] 


In  these  equations,  terms  with  combinations  of  c,  \j/,  V|/,2,  and  w,2  of  order 
5  and  higher  are  neglected  as  higher  order  terms . 

Because  in-plane  stress  and  strain  dominate  thin  beam  bending, 
through  the  thickness  shear  strain  (£4)  is  developed  using  only  the 
linear  terms  of  Eqn(2-22) .  Combining  these  linear  terms  with  Egn(20) 
and  Eqn(23)  yields  the  following  expression  for  through  the  thickness 
shear  strain: 


„  if  „  Zn  «2 

84  —  2823  —  ~  M3  2  i"(l  )^2,3 

i_£L  r  r] 


(2-21) 


Small  angle  approximations  are  used  in  the  global  displacement  equations 
because  shear  strain  is  linear,  resulting  in: 

M,  =  0 

f  -7^ 

(2-28) 


M2(j,z)  =  v|^1— ^j  +  z(\|/)+z^^(V|/  +w,2 ) 


M3 (j) =  w 

Substituting  Eqn(2-28)  into  Eqn(2-27)  using  k  =  4/(3h^),  one  obtains: 

1 


1-^ 


4z  4  3  8z  / 

^.2  -- :r-(v  +  w2)+:n:^(¥  +  W2) 

h  3h  r 


(2-29) 
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The  scale  factor  term  in  front  of  Eqn(2-29)  is  approximated  with  a 
truncated  binomial  expansion: 


1 


1- 


z 


=  1+^+ 
r 


r 


+. 


1 


The  is  also  assumed  negligible  since  z  <  h/2,  leaving: 

£4 


(2-30) 


where 


(2-31) 


£4  ='^.2+¥  and  X42 

as  the  expressions  for  through  the  thickness  shear  strain. 


2 . 5  Beam  Potential  Energy 

The  potential  energy  lip  is  defined  as  the  sum  of  the  internal 
strain  energy  (U)  and  the  external  work  (V) : 

np  =  f/+V  (2-32) 

The  internal  strain  energy  is  defined  as: 

U=jw’dV  (2-33) 

vol 

where  W*  is  the  strain  energy  density  function: 

!«(,«£«£«  (2-34) 


The  aijki  term  is  an  elasticity  tensor. 

In  this  case,  the  internal  strain  energy  is  made  up  of  two 
components,  the  in-plane  strain  (t/j)  and  the  through  the  thickness  shear 
strain  [1/2)  ■ 


k 

U,  =^bj  jQ^^eldzds 

^  I  -* 

2 


(2-35) 
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(2-36) 


ft 

U2  =\bj  ja^^eldzds 

^  I 

2 

Using  Eqn(2-25) : 

e^=(e“)^+2Xz'’X2pe2 +X(^'3C2.)'  (2-37) 

p=i  p=i 

Now  f/j  is  divided  into  three  parts,  each  representing  a  strain  component 
from  Eqn (2-37 ) : 

+^2+1^3)  (2-38) 

2  / 

where : 

k 

1^1  =  jQiki^if  ^z  =  A{t°f 

k 

"2 

h 

1^2  ~  J  ^Qlk  (®2  5C2p2^^^  ~  ^£2  5%2I  "^^^22  ■*'^X23  ■^^X24  "^(^25  "*■  ^X  26 

h  p=l 

~2 
fi 

^^3  ~  j  ^2k  ^X2p2^ 

_6  Vp='  J 

2 

~  ^21  +2£x2iX22  ■*■^(2X21X23  ■*■  X 22) ■*"2G(x 21X24  +  X 22X23 )i'2'^^(X 21X25  +X22X24)  (2-39) 
■^■^^(X  21X26  "^X  22X25  ■*■  X  23  X  24  )  ■*■  "^(^X  21 X  27  ■*■2X22X26  ■*■2X23X25  ■*^X24) 

■*^2^(X  22X27  ■'■X  23X26  ■*^X  24X25)  ■*■^(2X23X27  ■*■2X24X26  ■*■  X  25  )^*^2^(X  24X27  ■*^X  25X25) 

■'■^(2  x  25  X  27  ■*■  X  26  )  ■*■  2‘^X  26  X  27  ■*■  ^X  27 

The  [ij's  of  Eqn(2-39)  have  been  integrated  through  the  thickness  such 
that : 
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strain  energy  due  to  through  the  thickness  shear  is  similarly 


derived.  In  this  case: 

U2  =]-b\\l^ds 

^  t 

where:  (2-41) 

h 

=  je44J(e^"+2e4%42z"+(X42)'z']* 

_h 

2 

^AS{iilf+lDSzlx,,  +  FS{x,,f 

Again,  the  elasticity  terms  (AS,DS,FS)  are  found  by  integrating  through 
the  thickness : 


[AS,  DS,  FS]  =  I [1,  z\z^]dz  ( 2  -42 ) 

h 

"2 

The  strain  expressions  of  Eqn(2-25)  and  Eqn(2-31)  are  divided  into 
linear  and  nonlinear  parts  to  facilitate  placing  them  into  the  preceding 
strain  energy  formulations.  This  is  done  by  first  defining  the 
displacement  gradient  vector  d  such  that: 


d'^  =  {  V  V  2  w  w  2  W22  V  V  ,2 } 

Now  the  in  plane  strain  can  be  expressed  as: 


e“=L;d+idXd 

X2p=L>+|dXd 


(2-43) 


(2-44) 


where  p  =  1, 2, 3, 4, 5, 6, 7 .  The  Lj's  are  column  vectors  and  the  Hj-'s  are 
symmetric  matrices.  Combining  the  Lj's  and  Hj's  with  the  displacement 
gradient  vector  creates  the  linear  and  nonlinear  strain  terms, 
respectively.  This  is  demonstrated  for  the  mid  plane  strain  term  (62): 
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V 


''.2 

w 


e“={0  1  -c 


0 


w. 


0  0  0  W 


”'.2 

^.22 

V 

V,2 


0 

0 

c 

0 

0 

o' 

■  V  ■ 

0 

1 

—c 

0 

0 

0 
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—c 
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0 

0 

0 

w 

c 

0 

0 

1 

0 

0 

0 

^,2 

0 

0 

0 

0 

0 

0 

0 

W,22 

0 

0 

0 

0 

0 

0 

0 

V 

0 

0 

0 

0 

0 

0 

0 

y.2. 

(2-45) 


The  through  the  thickness  shear  strain  (£4)  expressions  are  less 
complicated  due  to  the  previous  assumptions  about  linearity  and  small 
angle  approximations : 


e?=s;d 

X42=S[d 


(2-46) 


Appendix  A  contains  the  complete  H,  L,  and  S  matrices. 

Miller's  derivation  of  the  Hj-  matrices  contains  non  linear  terms 
of  \|/  whereas  Creaghan's  Hj's  contained  only  constants.  Miller  treats 
the  Hj- ' s  as  constants  to  maintain  similarity  with  Creaghan  and  Palazotto 
and  Dennis.  This  linearization  is  accomplished  by  substituting  the 
previous  incremental  value  of  \|f  into  the  current  Hj-'s  (3;  13;  19). 
Although  this  assumption  simplifies  the  impending  variational  calculus, 
it  limits  the  accuracy  of  Miller's  theory  because  the  behavior  of  the 
bending  angle  is  not  properly  modeled  and  leads  to  inexact  equilibrium 
equations.  This  may  be  a  major  reason  the  theory  is  limited  to  moderate 
rotations  (46  degrees) .  Subsequent  research  should  investigate  this 
assumption. 


2-18 


The  total  beam  strain  energy  can  now  be  expressed  as: 


^  I 

where  K  is  a  matrix  of  constants,  N,  is  a  matrix  of  the  linear 

functions  of  d,  and  Nj  is  a  matrix  containing  quadratic  terms  from  d. 

This  expression  is  formed  by  substituting  Egns(2-44)  and  (2-46)  into 
Eqns(2-38),  (2-39),  (2-41),  and  (2-42)  as  outlined  by  Palazotto  and 

Dennis  (18) .  This  form  of  the  equation  simplifies  finite  element 
formulation.  As  previously  stated,  the  higher  order  terms  of  \|i  are 
treated  as  constants  in  Nj  and  Nj .  The  K  ,  Nj  ,  and  matrices  are 
listed  in  Appendix  A. 

2  ■  6  Finite  Element  Foirmulation 

The  beam  finite  element  used  for  the  current  formulation  is  shown 
in  figure  2-9.  The  element  has  three  nodes  with  degrees  of  freedom 
(DOF)  V,  w,  W2,  V,  and  1(^2  at  the  end  nodes  and  only  v  at  the  middle 
node.  This  element  is  similar  to  Miller's  except  for  the  addition  of 
the  1(^2  degree  of  freedom  at  the  end  nodes. 

The  additional  DOF  facilitates  Hermitian  shape  functions  for  the 
bending  angle  (\(/)  .  This  element  maintains  continuity  for  w  and  \(i 
and  their  derivatives  and  C®  continuity  for  v.  Previous  efforts  only 
maintained  C®  continuity  for  \|/. 

Adding  two  DOF  also  mandates  additional  boundary  conditions.  This 
element  experiences  clamped,  hinged,  free,  and  symmetric  boundary 
conditions  and  \|/2  must  be  defined  for  each.  Because  a  physical 


N,(d)  ^  N,(d) 


d  ds 


(2-47) 
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interpretation  of  Vj/j  is  impossible,  the  values  for  \j/2  under  these 
conditions  are  derived  mathematically  from: 

d^w  M 


"  ds^  “  El 


(2-48) 


Based  on  this  assxamption,  \(f2  =  0  at  any  boundary  condition  incapable  of 
supporting  a  bending  moment.  Table  2-2  lists  the  resulting  boundary 
conditions  used  for  the  DOF  at  nodes  1  and  2 .  For  the  current  work,  v 
at  node  3  is  never  specified. 


Table  2-2  Specified  Boundary  Conditions 


V 

w 

W,2 

¥ 

¥.2 

clamped 

0 

0 

0 

0 

free 

hinged 

0 

0 

free 

free 

0 

free 

free 

free 

free 

free 

0 

point  of  symmetry 

0 

0 

0 

0 

free 
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The  values  of  the  DOF  q(ri)  at  any  point  can  be  expressed  in 


natural  coordinates  in  terms  of  the  element's  DOF  using: 


V 

■  G, 
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Ga 
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Hu 
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0 
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0 

0 

0 

Hu., 
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0 

0 

0 

0 

Hn., 

H22.,_ 

i|/« 


(2) 


w 


w 


(2) 


where  the  shape  functions  are: 

1 


(2,  =-(11  -n)  G2=-(ii  +11)  G3  =  (i-ii  ) 

=^(2-3ti+ti^)  //,2  =^(1-T1-T1^  +r|^) 
H,,  =^(2  +  3ii-ti^)  =£(-1-ii+ii1  +ti1) 


(2-49) 


(2-50) 


The  displacement  gradient  vector  d  of  Eqn(2-43)  must  be  related  to  the 


DOF  before  using  the  beam  energy  equation,  Eqn(2-47) .  The  conversion  is 
done  in  natural  coordinates  using  the  following  equation: 


V|/« 
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(2-51) 
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where  D  is  the  shape  function  matrix  used  to  relate  the  d  to  the  DOF. 

By  multiplying  by  the  inverse  of  the  Jacobian  matrix  (J)  ,  Eqn(2-51)  is 
converted  to  global  coordinates. 

d(s)  =  J-'d(h)  =  J-'Dq  =  Tq  (2-52) 

where  J'*  is  defined  as : 


J-‘  = 


0 

J_ 

a 

0 

0 
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0 
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0 

0 
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J_ 
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where  a  is  half  the  element  length. 

Now  the  beam  strain  energy  from  Eqn{2-47)  based  on  the  DOF  is: 


-  N,  N, 

t/  =  -fefq^T'^ 

K+^+— ^ 

2  J 
^  1 

3  6 

(2-54) 


Integrating  along  the  beam  length  yields  new  terms  used  for  a  tangent 
stiffness  matrix: 


K  =  bjT'^kTds-Ni  =fejT^N,Trf5-N2  =  (2-55) 

/  I  I 

Substituting  Eqn(2-55)  into  Egn(2-54),  one  obtains  this  element's 
potential  energy  equation: 


N.(q) 

6 


q-q"R 


(2-56) 


where  R  is  a  vector  representation  of  external  forces . 

Defining  the  beam  equilibrium  equation  requires  taking  the  first 
variation  of  Eqn(2-56)  and  setting  it  equal  to  zero: 
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(2-57) 


5n^=5q^  K  +  q-R  =5q^[F(q)]  =  0 

where  F  is  the  resultant  equilibrium  equations.  The  first  variation  is 
accomplished  with  respect  to  the  vertical  displacement  (w)  because  the 
bending  angle  \|/  is  taken  to  be  a  linearized  function  and  the  membrane 
displacement  (v)  is  a  linear  function.  This  is  where  the  theory  goes 
wrong.  The  first  variation  should  be  accomplished  with  respect  to  the 
bending  angle  which  includes  a  third  order  approximation  as  well  as  the 
vertical  displacement.  The  linearization  of  \|/  limits  the  exactness  of 
the  equilibrium  equations  as  discussed  in  section  2.5. 

Because  5q  is  arbitrary  and  independent,  F(q)  =  0  becomes  the  set 
of  algebraic  nonlinear  equilibrium  equations.  To  solve  these  equations, 
Eqn(2-57)  is  expanded  in  a  truncated  Taylor  series 
such  that: 

F(q+8q)  =  F(q)+|^8q  =  0 

aq 

or:  (2-58) 

|^5q  =  -F(q) 

aq 

Substituting  the  expression  for  F  and  its  derivatives  into  Eqn(2-58) 
yields : 

N  N  ' 

KT8q  =  -  K+^+^  q  +  R 

where:  (2-59) 

Kt  =[K+N, +Nj] 

Eqn(2-59)  is  the  final  equation  used  to  find  the  equilibrium  values  for 
the  11  DOF.  In  any  case,  the  nonlinear  stiffness  matrices  vary  with 
load  and  displacement  throughout  the  structure. 
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2.7  Niimerical  Solution  Methods 


Modifying  Eqn(2-59)  to  include  all  elements  of  a  finite  element 
model ,  one  obtains : 

1  nf.  f .  iCtA  1 

q  +  R  (2-60) 

where  n  is  the  number  of  elements  in  the  model  and  j  is  the  integration 
over  the  length  of  any  particular  element .  The  global  load  vector  R 


n 

z 

ftjT''(K  +  Ni+N2) 

Tds 

8q  =  -I 

- 1 

H 

H 

Z'  ^  ^ 

Tds 

;=i 

1 

- 

J 

L  '  ' 

1  2  3  J 

contains  as  many  rows  as  DOF.  The  q  and  8q  vectors  are  global 
displacement  arrays  constructed  from  the  elemental  displacement  vectors. 

Eqn(2-60)  is  integrated  using  Gauss  quadrature.  Although  exact 
integration  would  require  7th  order  Gauss  quadrature,  5th  order 
quadrature  is  used  with  negligible  difference  (8;  13). 

Miller  employed  the  Newton-Raphson  solution  methods  which  are 
retained  for  the  current  work.  Either  the  displacement  control  method 
developed  by  Dennis  and  Palazotto  or  a  modified  Riks-Wempner  technique 
is  used  to  generate  equilibrium  paths  (13;  19;  27) .  Both  techniques  are 
briefly  reviewed  here. 

The  displacement  control  method  iterates  load  values  based  on 
incremental  displacements  until  equilibrium  is  achieved.  The  initial 
displacement  as  well  as  each  displacment  increment  are  prescribed.  This 
method  transverses  horizontal  limit  points  similar  to  point  A  in  figure 
2-10;  however,  vertical  limits  such  as  point  B  cannot  be  evaluated 
because  no  unique  solution  exists  at  that  displacement.  Figure  2-11 
graphically  depicts  the  displacement  control  method.  Starting  at  an 
equilibrium  point,  the  displacement  is  incremented  and  a  load  determined 
using : 


[K  +  N,(q..,)+N,(q:.,)]5q..,=- 


,  N,(q..,)  ^  N.(qt,) 


q^-i  +Ri 


(2-61) 
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Load 


Displacement 

Figure  2-10  Typical  Load  vs.  Displacement  Curve 


Load 


qr 

Displacement 

Figure  2-11  Displacement  Control 

Eqn{2-61)  comes  from  Palazotto  and  Dennis  (19).  Iteration  of  Eqn(2-61) 
is  continued  until  convergence,  for  a  given  displacement  increment, 
satisfies  equilibrium 

The  Tsai  and  Palazotto  modified  Riks-Wempner  technique  transverses 
limit  points  the  displacement  control  method  cannot  (27) .  This 
technique,  which  is  referred  to  as  the  Riks  method,  modifies  the 
equilibriiim  equation  such  that: 

K+— +— |q-A,R  (2-62) 

2  3  ) 
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where  X  is  a  user  specified  loading  parameter.  Eqn(2-62)  is  expanded 
for  a  solution  similar  to  Eqn(2-56),  resulting  in: 

K^8q,=5A,R-F(q,,X,)  (2-63) 

A  search  radius,  Al  is  defined  as  the  additional  constraint  required  for 
Eqn(2-63)  to  be  a  solvable  system.  This  value  is  defined  as: 

=  AqJ,,Aq,„  +AX^.„R"R  (2-64) 

The  search  radius  establishes  the  distance  of  the  new  equilibrium  point 
from  the  current  equilibrium  point.  Assuming  Al  is  small  enough  to 
accurately  capture  the  equilibrium  path,  it  may  be  approximated  using: 

A/^=Aq7,,Aq,,,  (2-65) 

This  definition  causes  an  unsymmtrical  global  stiffness  matrix. 

Reducing  Sq,  to  two  teirms  restores  the  required  symmetry: 

8q,  =8q„  +5X,8q,,, 

where  (2-66) 

5q,  =-K;'F(q„X,).  8q,,=-K;‘R 

Each  incremental  displacement  and  load  increment  values  are  based  on  the 

results  of  the  previous  increment  and  a  small  change  5: 

Aq,.^,  =Aq,+8q,^, 

AX(^i  =  AXj  +5Xj 

Substituting  Eqn(2-67)  into  Eqn(2-65)  yields: 

a8X]  +  bSk ,  +  c  =  0 

where 

a=5ql5q,^,  fe  =  25q^j(Aq, +5q,) 
c  =  (Aq,.  +  8q, )'"  (Aq,  +5q,.  )-Al^ 

The  preceding  system  of  equations  converges  when  Sq,-  becomes  smaller 


(2-67) 


(2-68) 
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than  a  user  defined  value.  After  convergence  at  load  step  m,  the  total 


displacement  q^,  the  total  load  parameter  "Kn,  and  the  incremental  load 
step  AZ„  are  calculated  using: 

=q™-i+Aq„ 

(2-69) 

m  m—i  m 

AL  =A/  , — ^ 

m  m— I  \j 


where  is  a  user  defined  number  of  iteration  prediction  and  N^.i  the 
number  of  iterations  required  for  convergence  of  step  m-1 .  Once  the 
search  radius  for  a  step  is  defined,  the  initial  load  parameter  is  found 
using: 


AX.|  =  + 


A/. 

(8qa8q,2)^' 


(2-70) 


2 ■ 8  Step  by  Step  Riks  Algorithm 

This  Riks  method  example  is  presented  by  Creaghan  and  Miller  and 
is  repeated  here  for  clarity  (3;  13).  Figure  2-12  graphically  depicts 
the  process.  For  each  Riks  solution,  the  estimated  nvimber  of  iterations 
{Ns) ,  total  number  of  steps  (m)  ,  the  incremental  displacement  limit 
(q,)  ,  and  the  initial  load  parameter  (Ao)  are  user  prescribed  parameters. 

1.  First  increment,  first  iteration:  compute  only  the  constant 
stiffness  matrix,  K. 

2.  First  iteration,  subsequent  increments:  compute  8q„=-K;'R. 

After  the  first  increment,  Kj  is  composed  from  q^.y  nonlinear 
displacement  terms . 

3.  Each  iteration,  each  increment:  compute  Aq,- =  AAi5q;2.  For  the  first 
increment,  AAi  =  Ao,  which  is  a  program  input.  Current  examples 
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Figure  2-12  Riks  Method 


4. 


5. 

6. 

7. 

8. 


9. 


used  A<)  =  0.1  as  did  Creaghan  and  Miller  with  acceptable  results  (3; 

13).  For  subsequent  increments,  AXi  =  Xi. 

Each  iteration,  each  increment:  compute  in  order: 

8q,=-K;'F(q,.?i,) 

5q,  =8q„ +8A,,8q,2  (2-71) 

Aq/+,  =  Aq,  +8q,^, 

Compute  A/ =  (AqJ^,Aq,+i)  . 

Update  with  q,- =  q^./ +  Aq,- . 

Solve  Eqn(2-68)  for  ±8A,;.  For  complex  roots,  return  to  step  2  and 
arbitrarily  adjust  A/^,  which  changes  AX] . 

Choose  ±8X,'  based  on  which  value  yields  a  positive  0  in: 

5q,  =8q,  ±8A,,8q,.2 

,  ,  (2-72) 

0  =(Aq,.+8q,.)Aq, 


If  both  are  positive,  then  bX,  =  -c/b. 

Update  the  displacement  and  loading  parameters  using: 

Aq,+,  =Aq,+8q.^, 

AX.(^j  =  AAij  +8A.j 


(2-73) 
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10.  Check  for  convergence.  If  no  convergence,  return  to  step  2.  upon 
convergence,  update  the  displacement,  load  parameter,  and  search 
radius  for  the  next  increment: 


=q„-i+Aq„ 

A/  —  A<  I  +  AA._ 

m  ffi— 1  m 


Al 


m 


=  Al 


m-1 


(2-74) 


where  is  a  user  defined  niomber  of  iteration  prediction  and  Nm.i 
the  n\amber  of  iterations  required  for  convergence  of  step  m-1. 

The  current  work  uses  Ns  =  2.5 . 

11.  Compute  the  loading  parameter  for  the  first  iteration  of  the  next 
increment : 

Al 

AX,  =± - 2^-;^  (2-75) 

(Sq^Sq^r 

12 .  Return  to  step  2  to  begin  the  next  increment . 


2 . 9  Arch  Geometry  Considerations 

As  previously  discussed,  geometry  plays  a  significant  role  in  the 
arch  behavior.  Figure  2-13  presents  a  typical  arch  geometry.  Miller 


S 


Figure  2-13  Arch  Geometry  Definitions 


considers  an  arch  thick  when  (12) : 
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(2-76) 


He  also  defines  a  deep  arch  using: 


A,>8 


where 


(2-77) 


These  classification  will  be  used  when  comparing  the  current  work  with 
Miller ' s . 
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III.  Results  and  Discussion 


3 ■ 1  Overview 

This  chapter  is  divided  into  two  sections.  In  the  first  section, 
the  current  work  is  compared  to  Miller's  work  to  assess  the  effects  of 
the  Hermitian  shape  functions.  The  second  section  presents  the  results 
of  extensive  arch  analysis  using  the  current  program. 

3 . 2  Comparison 

3.2.1  Clamped  Isotropic  Shallow  Thin  Arch.  The  first  problem 
considered  was  originally  presented  by  Belytschko  and  Glaum(l)  and  is 


Figure  3-1  Clamped  Isotropic  Shallow  Thin  Arch 


presented  here  in  figure  3-1.  Miller  and  Creaghan  also  analyze  this 
structure (3 ;  13).  Because  the  current  work  is  an  evolution  of  Miller's, 
comparisons  are  only  made  with  his  results. 

This  first  example  combines  large  displacements  and  geometric 
nonlinearity  with  small  rotations.  As  such,  it  is  used  to  validate  the 
current  work  and  demonstrate  that  the  updated  shape  functions  do  not 
corrupt  the  program's  capabilities.  The  structure  is  symmetrically 
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Figure  3-2  Comparison  for  Clamped  Shallow  Thin  Arch 


modeled  using  17  elements.  Miller's  model  contained  81  active  degrees 
of  freedom  (DOF)  while  the  current  work  used  99.  Both  formulations  use 
displacement  control  with  a  convergence  tolerance  of  0.1%  and  maximum  of 
200  iterations  per  displacement  increment.  Twenty  increments  are  used 
with  an  initial  displacement  of  2.54mm.  Convergence  is  based  on: 


^|IM: 

N 

xioo%<roL 


(3-1) 


where  q^rr  <X^r-ir  and  are  elements  of  the  global  displacement  vector, 
for  the  (r  -  and  the  first  load  increment,  i  =  1  to  n  degrees 

of  freedom,  and  TOL  is  the  user  defined  tolerance  (19) . 

Figure  3-2  compares  the  vertical  displacement  (w)  at  node  8  (the 
point  of  symmetry  where  P  is  applied)  for  each  solution  method.  The 
maximum  difference  is  less  than  1.0%  at  10.16mm  of  displacement.  This 
figure  is  attained  using: 


Miller' s- Current 
Miller' s 


difference  = 


xlOO 


(3-2) 


This  confirms  correct  introduction  of  Hermitian  shape  functions  into  the 
finite  element  formulation. 

Examining  the  derivative  of  the  bending  angle  (\|/,2)  demonstrates 
the  effect  of  the  higher  order  shape  function.  Figure  3-3  compares 
(v|/,2)  for  the  current  work  with  Miller's  results.  The  linear  shape 
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Figure  3-3  Clamped  Shallow  Thin  Arch  Rotation  Rate  of  Change  Comparison 


function  of  the  previous  effort  create  a  curve  with  step  discontinuities 
at  each  node.  The  current  work  provides  a  continuous  curve  for  the 
bending  angle  rate  of  change  expression.  Thus,  the  higher  order  shape 
function  is  indeed  properly  incorporated  and  does  not  degrade  the 
capabilities  of  the  previous  effort. 

To  examine  the  efficiency  of  the  two  programs,  the  number  of 
elements  used  to  model  the  structure  is  incrementally  reduced  until  the 
solutions  deviated  from  the  "approved"  solution  presented  in  figure  3-2 
by  more  than  1.0%.  Both  formulations  were  run  using  the  same  setup  as 
the  initial  solution  with  the  exception  of  the  reduced  niiinber  of 
elements.  Using  12  elements  and  56  active  DOF,  Miller's  program 
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deviated  0.848%.  Using  this  model,  his  program  required  23.1  seconds  to 
solve  using  a  Sun  Sparc20  workstation.  The  current  work  required  only  7 
elements  and  39  active  DOF  to  achieve  similar  results,  using  the  same 
convergence  parameters.  This  setup  required  only  17.3  seconds  to  run  on 
the  same  workstation,  a  25.1%  reduction. 

3.2.2  Cantilever  Isotropic  Thin  Beam.  Figure  3-4  presents  the 
second  problem  examined.  This  structure  experiences  large  deflections 
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Figure  3-4  Cantilever  Isotropic  Beam 


and  moderate  rotations.  Hsiao  and  Hou  initially  examine  this  beam  using 
an  updated  Lagrangian  formulation  (14) .  The  structure  is  modeled  with 
40  elements  providing  200  active  DOF  for  Miller's  solution  and  240  for 
the  current  program.  Displacement  control  was  used  with  an  initial 
displacement  of  0.15m,  58  increments,  a  limit  of  100  iterations  per 
increment,  and  a  tolerance  of  0.5%.  Figure  3-5  compares  the  tip 
displacement  (w)  solution  of  the  current  method  to  Miller's.  Because 
this  structure  experiences  large  rotations,  the  current  work  should 
deviate  from  Miller's  if  the  bending  angle  shape  functions  are  the 
limiting  factor  of  the  previous  work.  However,  the  two  solutions  are 
almost  identical.  This  is  the  first  indication  that  the  limit  of  the 
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Figure  3-5  Cantilever  Isotropic  Beam  End  Displacement  Comparison 


previous  effort  is  not  in  the  finite  element  formulation,  but  in  the 
derivation  of  the  equilibrium  equations. 

The  derivative  of  the  bending  angle  ('|/,2)  is  examined  for  the 
entire  beam  at  maximum  load.  Figure  3-6  presents  these  results.  As 
with  the  first  example,  the  current  effort  generates  a  continuous  curve 
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Figure  3-6  Cantilever  Beam  Rotation  Rate  of  Change  Comparison 


whereas  Miller's  results  are  not.  This  indicates  the  Hermitian  shape 
functions  do  properly  estimate  the  bending  angle  and  its  first 
derivative . 

The  efficiency  was  also  examined  using  the  procedure  outlined  in 
section  3.2.1.  Miller's  program  required  14  elements  and  70  active  DOF 
to  match  the  solution  generated  using  40  elements .  This  model  required 
45.4  seconds  to  solve.  The  current  work  required  fewer  elements  (13) 
but  more  DOF  (78)  .  This  caused  the  current  program  to  run  70.7  seconds 
before  reaching  a  solution. 

3.2.3  Very  Deep  Hinged  Clamped  Isotropic  Arch.  The  following 
example  combines  large  rotations  and  large  displacements.  The 
unsymmetrical  boundary  conditions  combine  with  the  very  deep  geometry, 
as  shown  in  figure  3-7,  to  present  a  formidable  test  of  any  large 


Figure  3-7  Very  Deep  Hinged  Clamped  Isotropic  Arch 


rotation  and  large  displacement  theory.  The  arch  was  modeled  with  100 
elements  and  the  displacement  control  method  was  used.  A  maximum  of  100 
iterations  was  specified  for  each  of  the  200  increments,  and  the 
tolerance  was  set  at  2.0%. 
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The  comparison  with  Miller's  work  is  shown  in  Figure  3-8.  The 
programs  produce  nearly  identical  results  up  to  the  collapse  load  of 


Figure  3-8  Load  vs.  Displacement  for  a  Very  Deep  Hinged  Clamped  Arch 


1.720  kn.  and  a  displacement  of  0.826  m.  At  this  point  the  solutions 
diverge.  However,  Miller  discovered  solutions  generated  after  collapse 
were  unreliable  due  to  large  energy  changes  within  the  elements  (13). 

As  such,  solutions  reached  after  collapse  are  neglected. 

The  rate  of  change  of  the  bending  rotation  (V|/,2)  is  presented  in 
figure  3-9  and  shows  the  same  results  as  the  previous  two  examples. 
Again,  continuity  of  the  bending  angle  is  attained.  Even  with  the 
higher  order  continuity,  the  current  model  does  not  provide  a  better 
solution.  These  results  clearly  indicate  a  limitation  caused  by  the 
linearization  of  the  rotation  angle  (i|/)  during  derivation  of  the  large 
rotation  kinematics . 

Although  the  rotation  limits  have  not  been  improved,  the  current 
work  does  address  another  shortcoming  of  the  previous  efforts.  At  large 
deflections.  Miller  noted  a  propensity  for  elements  near  the  hinged 
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Figure  3-9  Rotation  Rate  of  Change  Comparison 


boundary  to  "kink  over"  (13) .  Figure  3-10  shows  the  undeformed  and 
deformed  shapes  created  by  each  model  at  35  inches  of  vertical 
deflection  at  the  middle  node.  Despite  providing  similar  solutions  for 
the  displacement  at  the  middle  of  the  arch,  the  kinking  phenomenon 
Miller  noticed  is  nonexistent  in  the  current  solution.  As  a  result, 
displacements  near  the  hinged  boundary  are  significantly  different.  For 
this  reason,  the  current  model  should  be  used  for  investigating  large 
rotation  problems . 

Another  advantage  of  the  current  formulation  is  that  erroneous 
solutions  are  eliminated.  In  this  example.  Miller's  formulation 
continued  providing  solutions  after  the  model  experienced  element 
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Figure  3-10  Deformed  Shapes  for  Clamped  Hinged  Very  Deep  Arch 

kinking  despite  the  unstable  nature  of  the  structure  (13) .  The  current 
program  does  not  converge  on  a  solution  beyond  this  deflection  level. 

The  elimination  of  element  kinking  and  inaccurate  solutions 
clearly  demonstrate  the  advantage  of  using  Hermitian  shape  functions  for 
interpolation  of  the  bending  angle  (\|f)  when  modeling  arches  with  large 
displacements  and  large  rotations. 

For  this  example,  the  current  program  required  98  elements  with 
586  active  DOF  to  maintain  an  accurate  solution.  Miller's  program  only 
required  88  elements  and  440  active  DOF.  As  a  result,  the  current 
formulation  ran  for  1049  seconds  versus  805  seconds  for  the  previous 
effort.  Clearly,  the  trade  off  for  C^  continuity  is  a  finer  mesh  and 
larger  number  of  DOF,  which  in  turn  require  more  time  to  solve. 

3 . 3  Arch  Evaluation 

This  section  makes  use  of  the  newly  developed  curved  element  to 
characterize  certain  geometric  parameters.  Two  classes  of  composite 
arches  are  evaluated.  The  first  group  investigates  how  the  arc 
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length  {.S)  affects  the  collapse  load  (Pc)  .  The  second  group  evaluates 
the  effect  of  changing  the  thickness  ratio  {fi/r)  versus  Pc-  In  either 
case,  simply  supported  composite  arches  are  symmetrically  modeled  using 
the  current  program.  All  loads  were  applied  and  displacements  (w) 
measured  at  the  point  of  symmetry.  The  arches  are  constructed  of  a 
[0/90] s  lay-up  of  AS4-3501-6  graphite  epoxy.  Each  ply  is  0.25mm  thick, 
for  a  total  thickness  of  1mm.  The  lamina  material  properties  are: 

Ei=142  GPa;  E2=9 . 8  GPa;  Gi2=Gi3=6.0  GPa;  G23=4 . 8  GPa;  and  Vi2=0.3.  These 
arches  experience  large  deformations  and  moderate  rotations .  Riks 
method  is  used  for  all  simulations  with  a  tolerance  of  0.5%,  a  maximum 
of  500  increments,  and  a  limit  of  100  iterations  per  increment.  The 
initial  load  parameter  (^o)  was  set  at  0.1,  with  the  number  of 
iterations  estimated  as  2.5  {Ns),  and  a  maximum  allowable  load  change 
{K^)  of  2.0. 

3.3.1  Constant  h/r.  Variable  Arc  Length.  For  this  portion  of  the 
evaluation,  five  different  thickness  ratios  are  considered.  Each  h/r  was 
established  by  varying  the  radius  while  holding  the  thickness  at  1mm. 

For  each  h/r  the  arc  length  was  varied  to  determine  how  boundary  location 
affects  arch  behavior.  Table  3-1  presents  the  specifications  of  each 
arch  analyzed  for  this  section  along  with  the  depth  and  thickness 
parameters  used  for  evaluation.  Representative  figures  are  presented  in 
Appendix  B. 

Table  3-1  includes  two  classification  parameters  for  thickness  and 
depth  evaluation.  For  thickness,  either  h/r  or  c/h  may  be  used.  The 
c/h  parameter  takes  opening  angle  into  account  and  is  more  accurate  as 
the  opening  angle  increases.  Either  6/c  or  X  are  used  for  depth 

classification.  X  is  calculated  using  X  =  [12 (l-u^ )]• ^^ (r/h) • ^ (a/2 ) . 
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Table  3-1  Arch  Investigation  Matrix,  Constant  h/r,  h  =  0.001m 


h/r 

S  (m) 

r  (m) 

c/h 

8/c 

X 

0.1 

thick 

0.01 

0.01 

4.79 

thick 

0.26 

shallow 

2.94 

shallow 

0.1 

thick 

0.015 

0.01 

6.82 

thick 

0.39 

deep 

4.41 

shallow 

0.1 

thick 

0.025 

0.01 

9.49 

thick 

0.72 

deep 

7.36 

shallow 

0.1 

thick 

0.05 

0.01 

5.98 

thick 

3.01 

deep 

14.71 

deep 

0.05 

thick 

0.015 

0.02 

7.33 

thick 

0.18 

shallow 

3.12 

shallow 

0.05 

thick 

0.02 

0.02 

9.59 

thick 

0.26 

shallow 

4.16 

shallow 

0.05 

thick 

0.025 

0.02 

11.70 

thick 

0.32 

deep 

5.20 

shallow 

0.05 

thick 

0.03 

0.02 

13.63 

thick 

0.39 

deep 

6.24 

shallow 

0.05 

thick 

0.05 

0.02 

19.98 

thick 

0.72 

deep 

10.40 

deep 

0.05 

thick 

0.075 

0.02 

19.08 

thick 

1.36 

deep 

15.61 

deep 

0.05 

thick 

0.1 

0.02 

11.97 

thick 

3.01 

deep 

20.81 

deep 

0.01 

thin 

0.025 

0.1 

12.47 

thick 

0.06 

shallow 

2.33 

shallow 

0.01 

thin 

0.05 

0.1 

24.74 

thick 

0.13 

shallow 

4.65 

shallow 

0.01 

thin 

0.075 

0.1 

36.63 

thin 

0.19 

shallow 

6.98 

shallow 

0.01 

thin 

0.1 

0.1 

47.94 

thin 

0.26 

shallow 

9.31 

deep 

0.01 

thin 

0.15 

0.1 

68.16 

thin 

0.39 

deep 

13.96 

deep 

0.01 

thin 

0.2 

0.1 

84.15 

thin 

0.55 

deep 

18.61 

deep 

0.005 

thin 

0.05 

0.2 

24.93 

thick 

0.06 

shallow 

3.29 

shallow 

0.005 

thin 

0.075 

0.2 

37.28 

thin 

0.09 

shallow 

4.93 

shallow 

0.005 

thin 

0.1 

0.2 

49.48 

thin 

0.13 

shallow 

6.58 

shallow 

0.005 

thin 

0.15 

0.2 

73.25 

thin 

0.19 

shallow 

9.87 

deep 

0.005 

thin 

0.2 

0.2 

95.89 

thin 

0.26 

shallow 

13.16 

deep 

0.005 

thin 

0.3 

0.2 

136.3 

thin 

0.39 

deep 

19.74 

deep 

0.0025 

thin 

0.05 

0.4 

24.98 

thick 

0.03 

shallow 

2.33 

shallow 

0.0025 

thin 

0.075 

0.4 

37.45 

thin 

0.05 

shallow 

3.49 

shallow 

0.0025 

thin 

0.1 

0.4 

49.87 

thin 

0.06 

shallow 

4.65 

shallow 

0.0025 

thin 

0.15 

0.4 

74.56 

thin 

0.09 

shallow 

6.98 

shallow 

0.0025 

thin 

0.2 

0.4 

98.96 

thin 

0.13 

shallow 

9.31 

deep 

0.0025 

thin 

0.4 

0.4 

191.8 

thin 

0.26 

shallow 

18.61 

deep 

0.0025 

thin 

0.6 

0.4 

272.6 

thin 

0.29 

deep 

27.92 

deep 

Similar  to  the  thickness  classification  parameters,  as  a  increases,  the 
X  parameter  is  preferred  because  it  directly  accounts  for  the  opening 
angle . 

Figures  3-11  presents  a  typical  Pc  vs.  w  curve.  The  critical  load 
is  defined  as  the  local  peak  indicated  by  point  A.  In  the  cases  where 
no  local  peak  exists,  it  is  assiimed  the  structure  does  not  collapse  but 
has  post  buckling  features.  Figure  3-12  presents  the  Pc  v  S  data  for 
h/r  =  0.1.  As  arc  length  approaches  zero.  Pc  smoothly  approaches 
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Figure  3-12  Critical  Load  Comparison  for  h/r  =  0.1 


infinity  except  for  the  decrease  at  S'  =  0.01m.  Further  examination 
indicates  this  arch  experiences  local  snapping  and  does  not  totally 
collapse  (total  collapse  is  defined  as  global  snapping) .  That  is,  the 
arch  geometry  changes  but  is  immediately  able  to  support  a  load  with 
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minimal  deflection  increases.  Figure  3-13  presents  the  load 
displacement  curve  for  an  arch  with  arc  length  of  0.01m  and  thickness 
ratio  of  0.1.  As  indicated,  the  arch  initially  collapses  at  10.059  N  at 
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Figure  3-13  Load  v  Displacement  for  S  =  0.01m,  h/r  =  0.1 


initial  collapse  load  ^  10.56  N 


1 .7mm 


a  displacement  of  0.75inm.  This  arch  is  able  tosupport  the  same  load 
again  after  the  displacement  increases  only  1.7mm  to  2.45mm.  Figure  3- 
14  presents  the  actual  arch  shape  at  the  collapse  displacement.  The 
counterflexure  points  indicate  the  point  on  the  arch  where  the  radius 
changes  signs.  Outside  these  points,  the  arch  isconcave  downward  and 
between  the  counterflexure  points  the  arch  is  concave  upward.  The 
distance  between  the  points  is  indicative  of  the  stability  the  arch.  As 
the  distance  increases,  stability  decreases  along  with  the  ability  to 
support  a  load.  As  indicated  in  figure  3-14,  the  counterflexure  points 
of  the  S  =  0.01m,  h/r  =  0.1  arch  are  separated  by  1.5mm  or  15%  of  the  arc 
length.  This  small  separation,  combined  with  the  minimal  displacement 
change,  indicate  that  a  local  snapping  occurs. 
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Figure  3-14  Undeformed  and  Collapsed  Arch  Shapes;  S  =  0.01m,  h/r  =  0.1 

To  demonstrate  the  total  collapse  of  an  arch,  the  typical  load 
displacement  curve  for  an  arch  with  S  =  0.025m  and  h/r  =  0.1  is 
considered.  This  curve  is  presented  in  figure  3-15.  This  structure 
initially  collapses  at  a  load  of  10.76N  and  a  displacement  of  2.72mm. 


Figure  3-15  Load  v  Displacement  for  S  =  0.025m,  h/r  =  0.1 


3-14 


As  shown,  the  arch  does  not  recover  after  the  initial  collapse  point, 
indicating  a  global  snapping  or  total  collapse  of  the  structure. 

The  undeformed  shape  of  this  arch  is  presented  in  figure  3-16 
along  with  the  shape  at  collapse.  The  counter! lexure  points  are 


Figure  3-16  Undeformed  and  Collapsed  Arch  Shapes;  S  =  0.025m,  h/r  =  0.1 


separated  by  3.8mm  (15.2%  of  S)  at  collapse.  Both  this  arch  {S  =  0.025, 
h/r  =0.1)  and  the  previous  arch  {S  =  0.01,  h/r  =0.1)  have  similar 
counterflexure  separation  distances  (15%  of  S)and  collapse  loads 
(10. 5N) . 

Although  both  arches  collapse  at  similar  loads,  the  longer  arch 
displaces  almost  12%  of  its  length  before  collapsing  versus  7%  for  the 
shorter  arch.  Also,  the  longer  arch  does  not  recover  after  initial 
collapse,  indicating  a  that  global  snapping  occurs. 

Similar  results  were  observed  for  thiclcness  ratios  (h/r)  of  0.05, 
0.01,  0.005,  and  0.0025.  These  results  are  summarized  in  table  3-2  and 
presented  graphically  in  Appendix  B. 
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Table  3-2  Constant  h/r  Results 


h/r 

S  (m) 

Pc  (N) 

Wo  (mm) 

recover 

post  snap 
disp  (mm) 

type  of 
snapping 

0.1 

0.01 

10.587 

0.7471 

yes 

1.70292 

local 

0.1 

0.015 

12.616 

1.2727 

no 

n/a 

global 

0.1 

0.025 

10.761 

2.7196 

no 

n/a 

global 

0.1 

0.05 

6.2709 

7.059 

no 

n/a 

global 

0.05 

0.015 

4.8981 

0.8508 

yes 

1.9892 

local 

0.05 

0.02 

5.8055 

1.2548 

no 

n/a 

global 

0.05 

0.025 

5.6238 

1.7891 

no 

n/a 

global 

0.05 

0.03 

5.1597 

2.385 

no 

n/a 

global 

0.05 

0.05 

3.5331 

5.7071 

no 

n/a 

global 

0.05 

0.075 

2.7036 

12.055 

no 

n/a 

global 

0.05 

0.1 

1.7821 

14.187 

no 

n/a 

global 

0.01 

0.025 

0.3848 

0.66 

yes 

0.645 

local 

0.01 

0.05 

0.67085 

1.6038 

yes 

0.808 

local 

0.01 

0.075 

0.52102 

3.23 

no 

n/a 

global 

0.01 

0.1 

0.40145 

5.5168 

no 

n/a 

global 

0.01 

0.15 

0.26823 

11.795 

no 

n/a 

global 

0.01 

0.2 

0.19992 

19.769 

no 

n/a 

global 

0.005 

0.05 

0.20175 

0.97761 

yes 

2.31239 

local 

0.005 

0.075 

0.2405 

1.81 

yes 

5.6552 

local 

0.005 

0.1 

0.1984 

2.9864 

no 

n/a 

global 

0.005 

0.15 

0.13688 

6.4393 

no 

n/a 

global 

0.005 

0.2 

0.10279 

10.941 

no 

n/a 

global 

0.005 

0.3 

0.06824 

23.472 

no 

n/a 

global 

0.05 

0.04979 

yes 

0.6532 

local 

0.0025 

0.075 

0.07853 

yes 

2.8736 

local 

0.0025 

0.1 

0.08843 

1.654 

yes 

4.7291 

local 

0.0025 

0.15 

0.06757 

3.366 

no 

n/a 

global 

0.0025 

0.2 

0.05178 

5.829 

no 

n/a 

global 

0.0025 

0.4 

0.02605 

21.968 

no 

n/a 

global 

0.0025 

0.6 

0.0133 

47.795 

no 

n/a 

global 

The  results  presented  in  table  3-2  highlight  two  important  trends 
in  arch  behavior.  First,  the  arches  that  recover  after  collapse,  or 
experience  local  snapping,  are  generally  classified  as  shallow 
structures,  based  on  the  specifications  listed  in  table  3-1.  It  is 
interesting  to  note  that  for  smaller  radii  (r  =0.01m  and  0.02m),  the  8/c 
parameter  correctly  classifies  shallow  arches  as  those  that  experience 
local  snapping,  but  as  the  radius  increases,  the  A,  parameter  provides  a 
more  accurate  prediction. 
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The  second  trend  is  for  Pc  to  increase  as  h/r  increases  for  a 
given  arc  length.  This  indicates  the  influence  of  through  the  thickness 
shear,  as  will  be  described  in  the  next  section. 

3.3.2  Constant  Arc  Length,  Variable  h/r.  Table  3-3  presents  the 
specifications  of  each  arch  analyzed  for  this  section  along  with  the 
depth  and  thickness  parameters  used  for  evaluation. 


Table  3-3  Arch  Investigation  Matrix,  Constant  S,  h  =  0.001m 


S  (m) 

h/r 

r  (m) 

c/h 

8/c 

X 

0.025 

0.1 

thick 

0.01 

9.49 

thick 

0.72 

deep 

7.36 

shallow 

0.025 

0.05 

thick 

0.02 

11.70 

thick 

0.32 

deep 

5.20 

shallow 

0.025 

0.01 

thin 

0.1 

12.47 

thick 

0.06 

shallow 

2.33 

shallow 

0.05 

0.1 

thick 

0.01 

5.98 

thick 

3.01 

deep 

14.71 

deep 

0.05 

0.05 

thick 

0.02 

18.98 

thick 

0.72 

deep 

10.40 

deep 

0.05 

0.01 

thin 

0.1 

24.74 

thick 

0.13 

shallow 

4.65 

shallow 

0.05 

0.005 

thin 

0.2 

24.93 

thick 

0.06 

shallow 

3.29 

shallow 

0.075 

0.05 

thick 

0.02 

19.08 

thick 

1.36 

deep 

15.61 

deep 

0.075 

0.01 

thin 

0.1 

36.63 

thin 

0.19 

shallow 

6.98 

shallow 

0.075 

0.005 

thin 

0.2 

37.28 

thin 

0.09 

shallow 

4.93 

shallow 

0.075 

0.0025 

thin 

0.4 

37.45 

thin 

0.05 

shallow 

3.49 

shallow 

0.1 

0.05 

thick 

0.02 

11.97 

thick 

3.01 

deep 

20.81 

deep 

0.1 

0.01 

thin 

0.1 

47.94 

thin 

0.26 

shallow 

9.31 

deep 

0.1 

0.005 

thin 

0.2 

49.48 

thin 

0.13 

shallow 

6.58 

shallow 

0.1 

0.0025 

thin 

0.4 

49.87 

thin 

0.06 

shallow 

4.65 

shallow 

0.15 

0.01 

thin 

0.1 

68.16 

thin 

0.39 

deep 

13.96 

deep 

0.15 

0.005 

thin 

0.2 

73.25 

0.19 

shallow 

9.87 

deep 

0.15 

0.0025 

thin 

0.4 

74.56 

thin 

0.09 

shallow 

6.97 

shallow 

0.2 

0.01 

thin 

0.1 

84.15 

thin 

0.55 

deep 

18.61 

deep 

0.2 

0.005 

thin 

0.2 

95.89 

thin 

0.26 

shallow 

13.16 

deep 

0.2 

0.0025 
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Figure  3-17  presents  a  Pc  vs.  h/r  curve  for  arches  0.025m  long. 
This  curve  was  constructed  by  varying  the  radius  of  curvature  while 
maintaining  a  thickness  of  1mm.  As  shown  in  figure  3-18,  as  h/r 
decreases  due  to  the  increasing  radius,  the  depth  of  the  arch  decreases 
or  becomes  more  like  a  flat  beam. 
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generate  the  curve  are  considered  thick  based  on  c/h  <25.  As  such,  the 
deviations  observed  in  the  Pc  vs.  h/r  curve  demonstrate  the  nonlinear 
effects  caused  by  significant  through  the  thickness  shear  stresses 
experienced  by  thick  beams  and  arches. 

Figure  3-19  presents  the  Pc  vs.  h/r  for  an  arc  length  of  0.2m.  A 
slight  variation  in  the  critical  load  rate  of  change  was  observed.  In 
this  case,  all  three  S  =  0.2m  arches  were  thin  arches,  based  on  the 
c/h  >  25  criteria  presented  in  table  3-3.  As  such,  the  Pc  vs.  h/r  curve 


Figure  3-19  Pc  vs.  h/r  for  S  =  0.2m 


indicates  minimal  nonlinear  behavior  caused  by  through  the  thickness 
shear . 

The  critical  displacement  (wb)  for  these  arches  was  also  examined. 
Figure  3-20  shows  the  wb  vs.  h/r  curve  for  S  =  0.025m.  This  curve 
deviates  from  the  linear  path  as  indicated.  Figure  3-21  presents  the  wb 
vs.  h/r  curve  for  S  =  0.2m.  In  both  cases,  as  the  arches  modeled  become 
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deeper,  the  rate  of  change  of  the  critical  displacement  decreases, 
irregardless  of  the  presence  of  through  the  thickness  shear  effects. 

The  analysis  of  arches  with  constant  arc  length  and  thickness  and 
changing  radius  of  curvature  demonstrates  the  nonlinearities  introduced 
as  depth  and  through  the  thickness  shear  increase. 

Similar  vs.  h/r  and  wb  vs.  h/r  curves  were  constructed  for  arc 
lengths  of  0.05m,  0.075m,  0.1m,  and  0.15m.  These  curves  further 
exemplify  the  relationships  between  critical  load  and  thickness  and 
critical  displacement  and  depth.  Appendix  C  contains  these  curves. 
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IV.  Conclusions  and  Recoimnendations 


4 . 1  Summary  and  Conclusions 

This  effort  successfully  modified  the  existing  finite  element 
formulation  of  a  potential  energy  based  large  deformation  and  moderate 
rotation  theory.  Comparison  with  the  unmodified  formulation  indicated 
proper  incorporation  of  Hermitian  shape  functions  for  the  bending  angle 
(\|/)  .  Despite  the  higher  order  shape  functions,  the  rotation  limit  of 
the  current  work  remained  the  same  as  the  previous  effort.  Although  the 
limits  were  unchanged,  the  current  formulation  did  eliminate  element 
kinking  experienced  in  previous  efforts. 

Niimerous  arch  configurations  were  modeled  to  analyze  their 
behavior  with  respect  to  arc  length  and  thickness  variations.  As  the 
arc  length  decreased  for  a  given  thickness  ratio  (h/r) ,  the  critical 
load  (Pe)  approached  infinity  except  for  a  minor  decrease  as  the  arch 
became  flatter.  This  dip  in  the  Pc  vs.  S  curve  comes  about  because 
global  snapping  becomes  a  local  phenomena  before  snapping  disappears 
altogether.  The  decrease  in  critical  load  does  not  represent  a 
structural  failure,  but  a  geometry  change  which  does  not  affect  the  load 
carrying  capability  of  the  arch.  An  arch  experiencing  local  snapping 
recovers  its  structural  integrity  with  only  a  small  displacement 
increase.  On  the  other  hand,  after  global  snapping,  an  arch  experiences 
extremely  large  displacements  before  recovering.  The  extent  of  snapping 
(local  or  global)  can  be  accurately  predicted  by  examining  the  depth  of 
an  arch.  Deep  arches  experience  total  collapse  while  shallow  arches 
only  undergo  local  snapping. 
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Varying  the  radius  of  curvature  of  an  arch  while  holding  the 
thickness  and  arc  length  constant  demonstrated  independence  of  depth  and 
thickness .  Compared  to  an  analysis  where  through  the  thickness  shear  is 
neglected,  these  stresses  decreased  the  critical  load  of  both  shallow 
and  deep  arches  as  the  thickness  ratio  increased.  Similarly,  as  the 
depth  increased,  nonlinearities  decreased  critical  displacement  for 
thick  and  thin  beams  alike.  Thus,  the  critical  load  is  dependent  on  the 
thickness  and  the  critical  displacement  is  a  function  of  the  depth  of 
the  arch. 

4 . 2  Recommendations 

The  current  model  is  limited  to  45  degrees  of  bending  rotation 
just  as  the  previous  effort  was.  This  indicates  the  previous  finite 
element  formulation  was  sufficient  and  the  underlying  equilibrium 
equations  are  limiting  the  accuracy  of  the  program  beyond  the  bending 
rotation  limits.  As  such,  the  equilibrium  equations  should  be  rederived 
taking  the  first  variation  of  the  potential  energy  equation  with  respect 
to  the  bending  angle  (V|/)  as  well  as  vertical  displacement  (w)  .  Such  a 
derivation  could  then  be  used  with  the  current  finite  element 
formulation  to  increase  bending  angle  limit. 
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Appendix  A 

A.l  Qii  Transformations  of  Eqn(2-18) 

=  (2„  c»s'  e  +2((2,2  +2e^)sin^  0  cos"  sin'*  0 

Qn  =  {Qn  +Qn  -4^66  )sin'  9  cos"  0  +(2,2  (cos"  0  +sin"  0) 

Q22  =  <2,1  sin"  0  +  2((2,2  +  )sin"  0  cos"  0  +  cos"  0 
2,6  =  {Qn  -2,2  -2266)sin0  cos"  0  +(2,2  “222  +2266)sin'  9  cos9 
226  =  (2„  -2,2  -2266 )sin^  6  COS0  +(2,2-222  +2266)sin9  cos"  0 
266  =  {Qn  +222  -22,2  -2266 )sin"  6  cos'  6+266 (cos"  0  +  Sin"  0) 

244  =  255  sin"  0  +  244  cos"  0 

245  =(244  -255)cos0sin0 

2  55=  244  sin"  0  +  255  cos"  0 

A. 2  am  H  Matrices  in  Eqn(2-44)  and  Eqn(2-46) 
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A. 3  Element  Stiffness  Matrices  of  Eqn(2-47) 
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Appendix  B:  Critical  Load  Comparison  for  Constant  Thickness  Ratio  and 

Variable  Arc  Length 


Figure  B-1  Critical  Load  Comparison  for  h/r  =  0.05 


Figure  B~2  Critical  Load  Comparison  for  h/r  =  0.01 
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Appendix  C:  Critical  Load  and  Critical  Displacement  Comparisons  for 
Constant  Arc  Length  and  Variable  Thickness  Ratio 


Figure  C-1  Pc  vs.  h/r  for  S  =  0.05m 


Figure  C-2  Pc  vs.  h/r  for  S  =  0.075m 
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Figure  C-3  Pc  vs.  h/r  for  S  =  0.1m 


Figure  C-4  Pc  vs.  h/r  for  S  =  0.15m 
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Figure  C-5  vs.  h/r  for  S  =  0.05m 


Figure  C-6  vs.  h/r  for  S  =  0.075m 
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Appendix  D:  FORTAN  Program  Description 


D. 1  Overview 

The  following  program  was  initially  developed  by  Creaghan  and 
subsequently  modified  by  Miller  (3;  13).  The  attached  listing  does  not 
include  all  of  the  subroutines,  but  only  the  ones  modified  for  the 
current  work.  The  majority  of  the  modifications  were  done  to 
incoirporate  the  shape  function  changes  presented  in  Chapter  II.  Some 
minor  changes  were  also  incorporated  to  facilitate  multiple  simultaneous 
simulations . 

D.2  Data  Input  Format 

This  is  a  line  by  line  description  of  the  user  constructed, 
problem  specific  data  input  file  required  by  the  program.  Each  section 
represents  a  line  and  each  variable  on  a  line  is  separated  by  a  comma. 
Variables  are  identified  by  italics. 

1.  title:  problem  title  text  string 
2  .  linear,  isotro,  isarch,  i shape: 

a.  linear:  0  for  nonlinear,  1  for  linear 

b.  isotro:  0  for  composite,  1  for  isotropic 

c.  isarch:  0  for  straight  beam,  1  for  circular  arch 

d.  ishape:  0  if  half  a  symmetric  arch  is  modeled,  1  for  a  full 

arch  model 

3.  inctyp,  nic,  imax,  kupdate,  tol: 

a.  inctyp:  1  for  displacement  control,  2  for  Riks  method 

b.  nine:  number  of  increments  (load  or  displacement)  desired 

c .  imax:  maximum  n\amber  of  iterations  per  increment 
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4. 


5 


6. 

7. 

8. 

9. 


10. 


11. 

12. 


13. 


14. 

15. 


d.  kupdate:  not  used,  fill  with  0  or  1 

e .  tol :  convergence  tolerance  for  an  increment 

pincr,  eiter,  ttpi:  include  only  for  Riks  method  nonlinear  solutions 

a.  pincr:  initial  load  parameter,  Xo 

b.  eiter:  estimated  n^xmber  of  iterations  per  increment,  Ng  not  an 
integer 

c.  ttpi:  maximum  load  change  for  an  iteration  (X^ax) 

table {nine) :  include  only  for  nonlinear  displacement  control; 
this  is  a  list  of  nine  numbers  used  to  attain  the  incremental 
global  displacement 
nelem:  niJinber  of  elements 

delem{nelem) :  list  of  element  lengths;  1  for  each  element 
nbndry:  niimber  of  nodes  with  specified  boundary  condition  (BC) 
nboundi'nbndry,  5;  (v,V,\/, 2, w,w,2)  :  one  line  for  each  node  with  a 
specified  BC  the  first  number  is  the  node  number;  the  next  five 
numbers  describe  whether  or  not  the  dof  for  the  specified  node  are 
free  or  specified  in  the  order  indicated;  l=prescribed,  0=free. 
vbondiii) :  prescribed  displacements  for  the  fixed  dof  indicated  on 
the  previous  line 

Idtyp,  distld:  not  used;  fill  with  0,0.0 

nconc:  number  of  concentrated  loads  or  moments;  Riks  needs  at 
least  1. 

inconc {nconc)  :  use  only  if  nconoO;  dof  numbers  for  concentrated 
loads,  middle  node  dof  counts  here  (11  dof  per  element)  in  the 
following  order: 


v' , 'P' , w' ,  ,  v^  v" , 

vconc:  use  only  if  nconoO;  values  of  loads  specified  on  line  13 
ey,  nu,  ht,  width:  include  for  isotropic  material 
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a.  ey:  Young's  modulus 

b.  nu:  Poisson's  ratio 

c.  ht:  thickness 

d.  width:  beam  or  arch  width 

16.  el,  e2,  gl2,  nul2,  gl3 ,  g23,  width:  include  for  composite 
materials 

a.  el:  Young's  modulus  along  fibers 

b.  e2;  Young's  modulus  transverse  to  fibers 

c.  gl2:  shear  modulus  gi2 

d.  nul2:  Poisson's  ratio  V12 

e.  gl3:  shear  modulus  g^ 

f.  g23:  shear  modulus  g23 

g.  width:  beam  or  arch  width 

17.  nplies,  pthick:  include  for  composite  materials 

a.  nplies:  nximber  of  plies 

b.  pthick:  individual  ply  thickness;  same  for  all 

18.  theta(nplies) ;  include  for  composites;  list  of  ply  orientation 
angles  in  degrees 

19.  rad;  include  for  curved  arch;  arch  radius  of  curvature 

20.  nforc:  number  of  nodal  resultant  forces  to  be  calculated 

21.  i  force  (nforce)  :  include  only  for  nforoQ;  dof  numbers  for  force 
calculations,  include  middle  node  in  calculation 

22.  nstres:  not  used;  fill  with  0 
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D. 3  Program  Listing 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

0)  , 


program  beam 

See  bottom  of  file  for  variable  and  subroutine  listing. 

This  version  is  nl  n2  updated  to  include  tangent  function  and 
hermitian  shape  functions  for  the  rotation  d.o.f. 
remember  -e  compile  option  (will  give  132  character  lines) 
this  version  is  also  prompts  for  plot  and  shape  file  names 
implicit  double  precision  (a-h,o-z) 
character*64  gname, f name , pname , sname 


common/chac/gname, fname 

common/elas/ae,de, fe,he, ej , el, re, te, as, ds , fs 
common/input/tol , table (250) ,delem(250) , vbound(2500) ,distld, 

.  vconc ( 2  500 ) , ey , enu, ht , el , e2 , gl2 , enul2 , enu21 , gl3 , g23 , pthick, 
rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 

.  nelem, nbndry, nbound (250 , 6) , ldtyp,nconc, iconc (2500) , 

.  nplies,nforc, iforc (2500) ,nstres, istres (250) , ibndry (2500)  , 

.  theta(20) ,idload(250) ,coord(251) , width, nnod, pincr, eiter,  ttpi 

common/ St f/stif (11,11) , elp (11) ,eln(ll, 11) , eld (11) 

common/proc/gstif (2500,11) ,gn(2500,ll) ,gf(2500) ,gd(2500) , yperm  (250 
vpres(2500) 


call  r input (pname, sname) 
call  elast 

if (inctyp.eq.l)call  proces (pname, sname) 
if (inctyp.eq.2)call  rikspr (pname, sname) 
***********************************************************  ***** 


c 

C  VARIABLES  FOR  BSHELL 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


fname  input  file 
gname  output  file 
ae,de,  elasticity  terms 

fe,he,  elasticity  terms 

ej,el,  elasticity  terms 

re,te,  elasticity  terms 

as,ds,  elasticity  terms 

fs  elasticity  term 

ey  Young ' s  modulus  for  isotropic  case 

enu  Poisson's  ratio  for  isotropic  case 

ht  thickness  of  beam  for  isotropic  case 
el,e2,  laminate  material  properties 

gl2,enul2,  " 
enu21,gl3,  " 
g23 

pthick  laminate  ply  thicltness 

nplies  number  of  plies  in  laminate 

theta (20)  ply  orientation  angles 

tol  convergence  tolerance,  percent 

table (250)  displacement  increment  multiplicative 
factors 

delem(250)  element  lengths 

vbound(2500)  values  of  prescribed  displacement  boundary 
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c 
c 
c 
c 
c 
c 
c 
c 

increment 

c 


conditions 

distld  distributed  load  intensity 

vconc(2500)  concentrated  load  values 
rad  arch  radius  of  curvature 

linear  =1  for  linear  analysis,  =0  for  nonlinear 

isotro  =1  for  isotropic,  =0  for  laminate 

isarch  =1  for  arch,  =0  for  straight  beam 

ishape  =1  to  print  x,y  coordinates  for  each  node  at  each 


when  a  full  arch  is  represented  output  to  file  pname 
c  inclod  =1  to  increment  load(NA) ,  =0  increment  displacement 

c  nine  total  nximber  of  displacement  increments 

c  imax  maximum  number  of  iterations  per  increment 

c  nelem  total  number  of  elements  in  model 

c  nbndry  number  of  nodes  with  specified  boundary  conditions 

c  nbound(250, 6)  array  of  node  numbers  followed  by  I's  for 

c  fixed  b.c.'s,  zeros  for  unfixed 

c  Idtyp  =1  for  distributed  load,  =0  no  distributed  load 

c  nconc  total  number  of  concentrated  loads  input 
c  iconc(2500)  DOF's  for  specified  loads 

c  nforc  number  of  forces (including  moments) to  be  solved  for 
c  iforc(2500)  DOF's  at  which  to  calculate  forces 

c  nstres  number  of  elements  for  stress  calculation 

c  istres(250)  element  #'s  for  stress  calculation 

c  ibndry (2500)  DOF  numbers  for  b . c . ' s 

c  idload(250)  elements  with  distributed  load 

c  coord (251)  coordinate  of  the  nodes 

c  width  beam  or  arch  width 

c  nnod  number  of  nodes 

^*********************************************************************** 


c 

C  SUBROUTINES  FOR  BSHELL 

c 

c  rinput  reads  in  and  echos  input  data 

c  elast  computes  elasticity  terms 

c  proces  drives  the  solution  algorithm  for  displacement  control 

c  rikspr (pname, sname)  drives  the  solution  algorithm  for  Riks  method 

c  stiff  manages  stiffness  matrix  computations 

c  shape  computes  shape  function  array  dsf 

c  beamk  computes  constant  stiffness  array  bmk 

c  beamnl  computes  linear  stiffness  array  bmnl 

c  beamn2  computes  quadratic  stiffness  array  bmn2 

c  bndy  applies  displacement  boundary  conditions 

c  solve  solves  simultaneous  equations  in  banded  array  format 

c  converge  checks  solutions  for  convergence 

c  postpr  computes  nodal  loads  and  sends  to  output  file 

c 

end 

c 

c 

c 

subroutine  rinput (pname, sname) 
c 

character*  6  4  f name , gname 

character*4  title 

dimension  title (20) 

implicit  double  precision  (a-h,o-z) 
c 

common/chac/gname,  fname 
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common/ input /tol , table ( 250 ) , delem (250 ) , vbound ( 2500 ) , distld, 

.  vconc(2500) , ey , enu, ht , el , e2 , gl2 , enul2 , enu21 , gl3 , g23 ,pthick, 

rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 
nelem, nbndry,nbound (250 , 6) , ldtyp,nconc, iconc (2500 ) , 

.  nplies,nforc,iforc(2500) .nstres, istres (250) , ibndry (2500)  , 

.  theta (20) ,idload(250) , coord (251) ,  width, nnod, pincr, eiter, ttpi 

write (*, 1000) 

read (*, 1005) fname 

write (*, 1010) 

read ( * , 1005 ) gname 

write (*, 1011) 

read(*,  1005)pnaitie 

write ( * , 1012 ) 

read(*, 1005) sname 

open (5, file= fname) 

open ( 6 , f ile=gname, status= 'new' ) 

read(5,1015)title 

read ( 5 , * ) linear, isotro , isarch, ishape 
read  ( 5 ,  * )  inctyp ,  nine ,  imax,  Icupdte ,  tol 

i f ( 1 inear . eq . 0 . and . inctyp . eq . 2 ) read ( 5 , * ) pincr , eiter , ttpi 
if  ( linear. eq.  0  .and. inctyp. eq. Dread (5,  *)  (table  (i)  ,  i=l,ninc) 
read ( 5 , * ) nelem 

read (5,*) (delem(i) , i=l,nelem) 

calculate  nodal  coordinates 

nnod=nelem+l 
coord (1) =0 . 0 
do  5  ii=2,nnod 

coord(ii) =coord{ii-l) +delem(ii-l) 
read ( 5 , * ) nbndry 
do  10  i=l, nbndry 

read ( 5 , * )  (nbound ( i , j ) , j  =1 , 6 ) 
ifdof=0 

ifdof=counter  for  enumber  of  fixed  dof’s 

do  20  i=l, nbndry 
do  20  j=2,6 

if (nbound (i , j ). eq. 0) goto  20 
ifdof=ifdof+l 

ibndry (ifdof) = (nbound (i, 1) -1) *6  +  (j-1) 
continue 

read (5, *) (vbound (i) ,i=l, ifdof) 

read (5, *) Idtyp, distld 

i f ( Idtyp . eq . 1 ) read ( 5 , * ) ndload 

if (Idtyp. eq. 1) read (5 , *) (idload(i) ,i=l, ndload) 

read ( 5 , * ) nconc 

if (nconc .ne . 0) read(5 , * ) (iconc (i) , i=l, nconc) 

if (nconc .ne . 0) read (5 , *) (vconc (i) , i=l, nconc) 

if  ( isotro. eq. Dread (5,  * )ey, enu, ht, width 

if (isotro . eq. 0) read (5 , * ) el, e2 ,gl2 , enul2 , gl3 , g23 , width 

if (isotro . eq. 0) read(5 , *)nplies,pthick 

if (isotro . eq. 0) read (5 , * ) (theta (i) , i=l,nplies) 

if (isarch. eq. 1) read (5 , * ) rad 

read (5, *)nforc 

if (nforc .ne . 0) read (5 , * ) (if ore (i) , i=l,nforc) 
read(5, *) nstres 

if (nstres .ne . 0) read (5 , *) (istres (i) , i=l, nstres) 


c 

c  Echo  the  input  to  the  output  file 

c 

write (6,1015) title 

if (isarch.eq. 1) write (6, 1020) 
if (isarch.eq. 0) write (6, 1025) 
if (linear. eq.l)write(6, 1030) 
if (linear. eq.0)write(6, 1035) 
if (isotro . eq. 1) write (6,1040) 
if ( isotro. eq. 0) write (6,1045) 
if (ishape . eq. 1) write (6,1050) 
if (inctyp.eq. 1) write (6, 1060) 
if (inctyp.eq. 2) write (6,1055) 
write (6, 1065) nine 
write (6, 1070) imax 
write (6, 1075) tol 

if (inctyp.eq. 2) write (6, 1076) pincr,eiter,ttpi 
if (inctyp.eq. 1) write (6, 1078) 

if (inctyp. eq. 1) write (6,1080) (table (i) , i=l,ninc) 
write ( 6 , 1085 ) nelem 
write ( 6 , 1090 ) 

write (6, 10 95) (coord (i) ,i=l,nnod) 
write ( 6 , 1100 ) 
write (6, 1105) 
do  30  i=l,nbndry 

30  write (6, 1110) (nbound(i, j) , j=l,6) 
write ( 6 , 1115 ) ifdof 
write (6, 1120) (ibndry(i) ,i=l, ifdof) 
write (6, 1095) (vbound(i) , i=l, ifdof ) 
if (Idtyp.eq. 1) write (6, 1125)distld 
if (Idtyp . eq. 1) write ( 6 , 1130) (idload(i) , i=l, ndload) 
if (nconc .ne . 0 ) write ( 6 , 1135) 

if (nconc .ne . 0) write (6 , 1120) (iconc (i) , i=l, nconc) 
if (nconc .ne . 0 ) write (6,1095) (vconc (i) , i=l, nconc) 
if (isotro . eq. 1) write ( 6 , 1140) ey, enu,ht, width 
if (isotro . eq. 0) write (6 , 1145) el, e2 ,gl2 , enul2 , gl3 , g23 ,width 
if (isotro. eq. 0)write (6, 1150)nplies,pthick 
if (isotro. eq. 0)write (6, 1155) (theta(i) , i=l,nplies) 
if (isarch. eq. 1) write (6, 1160) rad 
write (6, 1165) (iforc (i) , i=l,nforc) 
write (6, 1170) (istres (i) , i=l,nstres) 
close (5) 
c  close (6) 

c 
c 

C  FORMATS 

c 

1000  f ormat (' Enter  your  input  file  name.') 

1005  format (A) 

1010  f ormat (' Enter  your  output  file  name.') 

1011  f ormat (' Enter  your  plot  file  name.') 

1012  f ormat (' Enter  your  shape  file  name.') 

1015  format (20a4) 

1020  format (/, lx, ' Element  type:  arch') 

1025  format (/, lx, ' Element  type:  straight  beam') 

1030  format (/, lx, 'Analysis  type:  linear') 

1035  format (/, lx, 'Analysis  type:  nonlinear') 

1040  format (/, lx, 'Material  type:  isotropic') 

1045  format (/, lx, 'Material  type:  laminate') 

1050  f ormat (/, lx, 'Printout  of  nodal  x,y  coordinates  requested') 
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1055  format (/, lx, 'Riks  method  specified') 

1060  format(/, lx, 'Displacement  control  method  specified') 

1065  format ( / , lx, ' Increments  specified: ' , 2x, i3 ) 

1070  format (/, lx, 'Maximum  iterations  specif ied: ', 2x, i3 ) 

1075  format (/, lx, ' Percent  convergence  tolerance :', 2x, dl2 . 5 ) 

1076  format (/, lx, 'pincr=' , 2x,dl2 .5,2x, •eiter=' , 2x, dl2 . 5 , 2x, 

'ttpi=' ,2x,dl2 .5) 

1078  format (/, lx, 'Displacement  Increment  Table') 

1080  format (8 (2x, dl2 . 5) ) 

1085  format  (/,  lx, 'N\Jimber  of  elements 2x,  i3 ) 

1090  format {/, lx, 'Nodal  Coordinates:') 

1095  format (8 {2x,dl2 . 5) ) 

1100  format (/, lx, 'DISPLACEMENT  BOUNDARY  CONDITIONS,  1=PRESCRIBED, 
X0=FREE' ) 

1105  format (/,4X, 'NODE  V  PSI-S  PSI-SS  W  W-S  ') 

1110  format (4x, i4 , 3x, 5 (i3 , 3x) ) 

1115  format (/, lx, 'NUMBER  OF  PRESCRIBED  DISPLACEMENTS:', 

.  i5, /, lx, 'SPECIFIED  DISPLACEMENT  DOF  AND  THIER 

.  VALUES  FOLLOW: ' ) 

1120  format (16i5) 

1125  format (/, lx, 'Distributed  Load  Intensity: ' ,2x,dl2. 5) 

1130  format(/, lx, 'Elements  with  distributed  load: ' , /, lx, 16i5) 

1135  format (/, lx, 'DOF  and  specified  concentrated  loadsfollow: ' ) 

1140  format!/, lx, 'Isotropic  material  properties  ey,  enu,  ht,  width:' 
, / , lx, 4dl2 .5) 

1145  format!/, lx, 'Composite  material  properties  el,  e2,  gl2,  enul2, 
gl3,g23,  width: ' , / , lx, 7dl2 . 5) 

1150  format !/, lx, 'Number  of  plies :', 2x, i3, 2x, ' Ply  thickness :', 2x, 

.  dl2.5) 

1155  format !/, lx, ' Ply  orientation  angles : ' , / , lx, 8 !2x, dl2 . 5) ) 

1160  format!/, lx, 'Radius  of  curvature: ', 2x, dl2 . 5) 

1165  format!/, lx, 'DOFs  for  equivalent  load  calculation:',/, 

,  lx,16i5) 

1170  format!/, lx, 'Elements  for  stress  calculation: ',/, lx, 16i5) 

1175  format !/, lx, i5) 
return 
end 
c 
c 
c 

subroutine  elast 
c 

implicit  double  precision  !a-h,o-z) 
c 

dimension  qbar !3 , 3 ) , rtheta !20) 
c 
c 

character*64  gname, fname 
common/ chac / gname , fname 

common/elas/ae, de, fe, he, ej ,el, re, te, as, ds , f s 

common/ input /tol , table !250 ) , delem!250 ) , vbound !2500 ) , distld, 

.  vconc!2500) ,ey,enu,ht,el,e2,gl2,enul2,enu21,gl3,g23,pthick, 

rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 
nelem,nbndry,nbound!250, 6) , ldtyp,nconc, iconc !2500) , 

.  nplies,nforc, iforc !2500) ,nstres, istres !250) , ibndry !2500)  , 

.  theta !20) ,idload!250) , coord! 2 51) ,width,nnod,pincr, eiter, ttpi 
c 
c 

c  Isotropic  case 
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c 

c  write { 6 , 1000 ) el , e2 , gl2 , enul2 , enu21 , gl3 , g23 , pthick 

if (isotro . eg. 0) goto  100 
gs=ey/ (2* (1+enu) ) 

denoin=l .  -enu**2 
qll=ey/denom 
ql2=enu*ey/denom 
q22=qll 

q2hat=q22- (ql2**2/qll) 

qs4=gs 

ae=q2hat*ht 

de=q2hat*ht**3/ (3*2 . **2) 
fe=q2hat*ht**5/ (5*2.**4) 
he=q2hat*ht**7/ (7*2 . **6) 
ej=q2hat*ht**9/ (9*2 . **8) 
el=q2hat*ht**ll/ (11*2 . **10) 
re=q2hat*ht**13/ (13*2.**12) 
te=q2hat*ht**15/ (15*2. **14) 
as=qs4*ht 

ds=qs4*ht**3/ (3*2 . **2) 
fs=qs4*ht**5/ (5*2 . **4) 
goto  200 
c 

c  Laminate  case 
c 

100  ht=pthick*nplies 

enu2 l=e2 *enul2 /el 
denom=l . -enul2  *enu21 
qll=el/denom 
ql2=enul2*e2/denom 
q22=e2 /denom 
c 

************************************************* 
c  calculate  the  elasticity  matrices  * 

c  * 

c  remem  that  the  z  axis  points  down,  * 

c  however,  the  first  ply  is  the  top  ply,  ie,  * 

c  the  ply  with  the  most  negative  z  ! 1 !  * 

Q**********  ************  **********************  ****** 

c 

c  initialize  elasticity  terms 

c 

ae=0 . 

de=0 . 
fe=0. 
he=0 . 
e  j  =0 . 
el=0. 
re=0 . 

te=0 . 
as=0 . 
ds=0 . 
fs=0 . 

do  45  ii=l,nplies 

45  rtheta (ii) =theta (ii) *3 . 14159265/180 . 

do  50  kk=l,nplies 

qbar (1, 1) =qll* (cos (rtheta (kk) ) **4) +2*ql2* (sin (rtheta (kk) ) **2) * 
(cos (rtheta (kk) ) **2 ) +q22* (sin(rtheta (kk) ) **4) 
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qbar (l,2)={qll+q22) * (sin (r theta (3ck) ) **2) * (cos (rtheta(kk) ) **2) + 
ql2* (sin(rtheta{kk) ) **4+cos (rtheta (kk) ) **4) 


qbar (2,2) =qll* (sin (rtheta (kk) ) **4)+2*ql2* (sin (rtheta (kk) ) **2) * 
(cos (rtheta (kk) )**2)+q22*cos (rtheta (kk) )**4 
qs4=gl3*dcos (rtheta (kk) ) **2+g23*dsin (rtheta (kk) ) **2 
q2hat=qbar ( 2 , 2 ) - ( qbar ( 1 , 2 ) *  *2 / qbar (1,1)) 
zl=(kk*l.  -  nplies* . 5) *pthick 
zu=zl-pthick 
ae=ae  +  q2hat*pthick 
de=de  +  q2hat* (zl**3-zu**3 ) /3 . 

fe=fe  +  q2hat* (zl**5-zu**5) /5 . 
he=he  +  q2hat* (zl**7-zu**7) /7. 
ej=ej  +  q2hat* (zl**9-zu**9) /9. 
el=el  +  q2hat* (zl**ll-zu**ll) /ll. 

re=re  +  q2hat* (zl**13-zu**13) /13 . 
te=te  +  q2hat* (zl**15-zu**15) /15 . 
as=as+qs4*pthick 

ds=ds+qs4* (zl**3-zu**3) /3 . 
fs=fs+qs4* (zl**5-zu**5) /5 . 

50  continue 

c  200  open(6, fiel=gname, status='old' ) 

200  write (6, 1000) ae, de, fe, he, ej ,el, re, te, as, ds, fs 
c  close(6) 

1000  format (/, lx, 'Elasticity  terms : ’ , / , lx, 8 (2x, dl2 . 5) ) 
return 
end 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

0)  , 
c 


1 


subroutine  proces (pname, sname) 
implicit  double  precision  (a-h,o-z) 
character*64  gname,fname 


common/chac/gname,  fname 

coramon/elas/ae, de, fe,he,ej ,el,re,te,as,ds,  fs 

common /input/ to 1, table (250) , delem(250) , vbound (2500) , distld, 
vconc(2500) ,ey,enu,ht,el,e2,gl2,enul2,enu21,gl3,g23,pthick, 
rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 
nelem,nbndry,nbound(250, 6) , ldtyp,nconc, iconc (2500) , 
nplies, nforc, iforc (2500) ,nstres, istres (250) , ibndry (2500)  , 
theta (20) ,idload(250) , coord (251) ,width,nnod,pincr,eiter, ttpi 
common/stf /stif (11,11) , elp (11) , eln(ll, 11) , eld (11) 

common/proc/gstif (2500,11) ,gn(2500,ll) ,gf(2500) ,gd(2500) , vperm(250 

vpres(2500) 

ndo  f  =nnod*  5  +ne 1 em 

ncount=l 

icount=l 

do  1  ii=l,ndof 

gd(ii) =0 . OdO 

do  2  ii=l,nbndry*6 
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vpres (ii) =0 . OdO 
vperm { i i ) =vbound ( i i ) 

start  new  increment  or  iteration/ 

zero  out  global  stiffness  matrices  and  global  force 
vector 

do  5  ii=l,ndof 
gf (ii) =0 . OdO 
do  5  jj=l,ll 
gstif (ii, j j)=0.0d0 
gn(ii, jj)=0.0d0 
kcall=0 

increment  prescribed  displacement  for  displacement  control 

if (linear. eq.l) goto  9 
if (icount .ne . 1) goto  9 
do  7  ii=l,nbndry*6 

if (ncount .eg. 1) vbound(ii) =vperm(ii) * table (1) 
if (ncount .gt . 1) vbound(ii) =vperm(ii) * (table (ncount) - 
table (ncount-1) ) 

loop  over  all  elements  for  stiffness  and  forces 

do  30  ielem=l,nelem 

do  10  ii=l,ll 

eld(ii) =gd(ii+ (ielem-1) *6) 

lccall=kcall+l 

call  stiff (ielem, icount, ncount, kcall) 

Assemble  global  stiffness  array,  gstif,  global  equilibrium 
stiffness,  gn,  in  banded  form.  Half-bandwidth=ll .  Also 
assemble  global  force  vector,  gf. 

nr= (ielem-1) *6  +  1 
do  30  jj=0,10 

gf (nr+ j j ) =gf (nr+ j j ) +elp ( j j +1 ) 
do  30  kk=l,ll-jj 

gstif (nr+j j ,kk) =gstif (nr+j j ,kk) +stif ( j j+l,kk+j j ) 
if (linear . eq. 1) goto  30 

if (icount. eq.l  .and.  ncount. eq.l) goto  30 
gn(nr+j j ,kk) =gn(nr+j j ,kk) +eln ( j j+1, kk+j j ) 
continue 

impose  force  boundary  conditions 
at  this  point,  gf=R 

if (nconc.eq.O)goto  45 
do  40  ii=l,nconc 
nb=iconc (ii) 
gf (nb) =gf (nb) +vconc (ii) 
continue 

calculate  the  residual  force  vector  for  nonlinear 
analysis .  - [gn] * {gd}+R=- [k+nl/2+n2/3 ] * {q}+R=gf 

if (icount. eq.l) goto  65 
do  60  ii=l,ndof 


add=0 . 

do  50  kk=l,ii-l 
if(ii-kk+l  .gt.  ll)goto  50 
add=add+gn(kk, ii-kk+1) *gd(kk) 

50  continue 

res=0 . 

do  55  jj=l,ll 

if(jj+ii-l  .gt.  ndof)goto  55 
res=res  +  gn (ii , j j ) *gd( j j+ii-1) 

55  continue 

c 

c  add  to  existing  gf  which  already  contains  R 

c 

gf (ii) =gf (ii) -res-add 
60  continue 
65  continue 
c 

c  impose  displacement  boundary  conditions 

c 

i f  ( icount .  eq .  1 )  cal  1  bndy  (ndof ,  gsti f ,  gf ,  nbndiry ,  ibndry ,  vbound) 
if (icount.gt.l)call  bndy (ndof, gstif,gf,nbndry, ibndry, vpres) 
c 

c  solve  system  of  equations,  result  in  gf 

c 

call  solve (ndof , gstif, gf, 0 , detm, detml) 
c 

c  update  total  displacement  vector  gd 

c 

do  70  ii=l,ndof 
70  gd(ii) =gd(ii) +gf (ii) 
if  (linear. eq.Dgoto  80 

call  converge (ndof , neon, icount, tol, imax) 
c 

c  if  no  convergence  (ncon=0)  start  next  iteration 

c 

if (neon. eq. 0) goto  3 
80  continue 

if(ncon.eq.l  .and.  ncount.le. nine) then 

call  postpr (icount, ncount,kcall, ndof, pname, sname) 
if (ncount .eq. nine) stop 
ncount=nco\mt+l 
icount=l 
goto  3 
endif 
return 
end 
c 

subroutine  bndy (ndof , s , si , ndum, idum, vdum) 
c 

C  . .  •  •  ■ 

c  subroutine  used  to  impose  boundary  conditions  on  banded  equations 

c  . 

c 

implicit  double  precision  (a-h,o-z) 

dimension  s  (2500 , 11) , si (2500) 

dimension  idum(ndum*6) , vdum(ndum*6) 

do  300  nb  =  1,  ndum* 6 

ie  =  idum(nb) 

sval  =  vdum(nb) 

it=10 
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i=ie-ll 

do  100  ii=l,it 
i=i+l 

if  (i  .It.  1)  go  to  100 
j=ie-i+l 

si (i) =sl (i) -s (i, j ) *sval 
s (i, j)=0.0 
100  continue 

s (ie, 1) =1 . 0 
si (ie) =sval 
i=ie 

do  200  ii=2,ll 
i=i+l 

if  (i  .gt.  ndof)  go  to  200 
si (i) =sl (i) -s (ie, ii) *sval 
s (ie, ii) =0 . 0 
200  continue 
300  continue 
return 
end 
c 
c 
c 
c 

subroutine  converge  (ndof ,  neon,  icoxint,  tol ,  imax) 

. . .  ■  ■  ■ ; . 

c  checks  for  convergnece  using  global  displacement  criterion 

. . 

implicit  double  precision  (a-h,o-2) 

common/proc/gstif (2500, 11), gn (2500, 11), gf (2500), gd (2500) ,  vperm(250 

0)  , 

vpres(2500) 

c 

rcurr=0 . 
do  10  m=l,ndof 

10  rcurr=rcurr  +  gd(m)*gd(m) 
if (icount .eq. 1) rinit=rcurr 
if ( icount. eq.l)ncon=0 
if (icount. eq.l) goto  20 
c  new  criteria 

ratio=100.  *  abs (sqrt (rcurr) -sqrt (pvalue) ) /sqrt (rinit) 
if ( ratio. le. tol )ncon=l 
20  pvalue=rcurr 

write (*, 100) neon, ratio, rinit, rcurr 
100  format (lx, 'ncon=  ' , i3,3x, 'ratio=  ',dl4.6,'  rinit=  ',dl4.6, 

X  '  rcurr=  ',dl4.6) 

if (icount . eq. imax) write ( 6 , 200) 
if (icount . eq. imax) stop 
200  format (lx, ' icount  equals  imax') 
if (ncon.eq. 0) icount=icount+l 
return 
end 
c 
c 
c 

subroutine  rikspr (pname , sname) 
c 

implicit  double  precision  (a-h,o-z) 
c 

character* 64  gname,fname 
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c 

c 

c 


common/chac/gname, fname 

common/elas/ae,de, fe,he, ej ,el, re, te, as, ds, f s 


common/input/tol, table (250) ,delem(250) ,vbound(2500) ,distld, 

.  vconc(2500) ,ey,enu,ht,el,e2,gl2,enul2,enu21,gl3,g23,pthick, 

rad, linear, isotro, isarch, ishape, inctyp, nine, imax, 

.  nelem,nbndry,nbound{250, 6) , Idtyp.nconc, iconc (2500 ) , 

nplies,nforc, iforc (2500) ,nstres, istres (250) , ibndry (2500)  , 

.  theta(20) , idload(250) ,coord(251) , width, nnod, pincr, eiter, ttpi 
common/stf /stif (11,11) , elp (11) , eln(ll, 11) , eld (11) 


c 

0)  , 


common/proc/gstif  (2500,11)  ,gn(2500, 11)  ,gf  (2500)  ,gd(2500)  ,.vperm(250 
vpres (2500) 


c 

dimension 

gld(2500) ,gld0(2500) ,gldl(2500) ,gdis(2500) , gstiOO (2500 , 11) , 
gf0(2500) ,gd00(2500) 
ndo  f =nnod*  5 +ne 1 em 
ncount=l 
icount=l 
iicut=0 

do  1  ii=l,ndof 

1  gd(ii)=0.0d0 

do  2  ii=l,nbndry*6 
vpres (ii) =0 . OdO 

2  vperm(ii) =vbound(ii) 
c 

c  start  new  increment  or  iteration/ 

c  zero  out  global  stiffness  matrices  and  global  force 

c  vector 

c 

tpincr=0 . 0 

if (ncount.eq.l)goto  2993 


c 

c 

c 

3 

2993 


2992 


c 

c 

c 


4 


5 


c 

c 

c 

c 

c 

c 


start  new  increment 

if  (iicut .  eq.  0)  dss=dss*eiter/icotint 
icount=l 

do  2992  ii=l,ndof 
gldO (ii) =0 . OdO 
gdOO (ii) =gd(ii) 

start  new  iteration 

do  5  ii=l,ndof 
gfO (ii)=0.0d0 
do  5  jj=l,ll 
gstif (ii, j j ) =0 . OdO 
gn (ii , j  j ) =0 . OdO 
]ccall=0 

increment  prescribed  displacement  for  displacement  control 


if  (linear. eq.Dgoto  9 
if (icount.ne.l)goto  9 
do  7  ii=l,nbndry*6 
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c  if (ncount . eq . 1 . and . iicut . eq . 0 ) vbound ( i i ) =vpenn  { i i ) *  table ( 1 ) 

c  7  if (ncount . gt . 1 . or . iicut . gt . 0 ) vbound ( ii ) =vperm ( ii ) * ( table (ncount ) - 
c  .  table (ncount-1) ) 

c 

c  loop  over  all  elements  for  stiffness  and  forces 
c 

9  do  30  ielem=l,nelera 
do  10  ii=l,ll 

10  eld(ii) =gd(ii+ (ielem-1) *6) 
c 

kcall=kcall+l 

call  stiff (ielem, icount, ncount, kcall) 
c 

c  Assemble  global  stiffness  array,  gstif,  global  equilibrium 
c  stiffness,  gn,  in  banded  form.  Half-bandwidth=ll .  Also 

c  assemble  global  force  vector,  gf. 

c  * 

nr= (ielem-1) *6  +  1 
do  30  jj=0,10 

gfO (nr+j j ) =gf0 (nr+j j ) +elp ( j j+1) 
do  30  kk=l,ll-jj 

gstif (nr+j j ,kk) =gstif (nr+j j ,kk) +stif ( j j+l,kk+j j ) 
if  (linear. eq.Dgoto  30 

if (icount .eq. 1  .and.  ncount. eq.l  .and. iicut. eq.O) goto  30 
gn (nr+j j , kk) =gn  (nr+j  j , kk) +eln ( j j+1, kk+j j ) 

30  continue 


c 

c  impose  force  boundary  conditions 

c  at  this  point,  gf=R 

c  /' 

if (nconc.eq.O)goto  45 
do  40  ii=l,nconc 
nb=iconc (ii) 

40  gf 0 (nb) =gf 0 (nb) +vconc (ii) 

45  continue 

do  47  ii=l,ndof 

47  gf (ii) =gf0 (ii) 
do  48  ii=l,ndof 
do  48  jj=l,ll 

48  gstiOO (ii, j j ) =gstif (ii, j j ) 

call  bndy (ndof,  gstif ,gf,nbndry, ibndry, vbound) 
c 

call  solve (ndof , gstif , gf, 0 , detm, detml) 

dss0=0 . 0 

do  49  ii=l,ndof 

gdis (ii) =gf (ii) 

49  dss0=dss0+gf (ii) *gf (ii) 
if (icount .ne . 1)  go  to  144 

detml=detm2 

detm2=detm 

if (ncount. eq.l. and. detm. It. 0.0  .and. iicut . eq. 0)  pincr=-pincr 
if (ncount. eq.l. and. iicut. eq.O)  dss=pincr*dsqrt (dssO) 
if (ncount. ne.l. or. iicut. gt.O)  pincr= 
dss/dsqrt (dssO) *detm*detml 

*pincrl/dabs (pincrl) 


c 

c  attempt  at  offloading  at  bifurcation  points 

c 

c  if (iicut. eq.l)pincr=-pincr 

c 
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pincrl=pincr 
prs=0 . 0 

do  142  ii=l,ndof 

142  prs=prs+gfO (ii) *gld(ii) 
stifpa=pincr*prs 

do  143  ii=l,ndof 

143  gld{ii) =pincr*gdis (ii) 

144  continue 
c 

c  calculate  the  residual  force  vector  for  nonlinear 

c  analysis.  - [gn] * {gd}+R=- [k+nl/2+n2/3] * {q}+R=gf 

c 

if (icount.eq.l)goto  69 
do  60  ii=l,ndof 
add=0 . 

do  50  kk=l,ii-l 
if(ii-kk+l  .gt.  11) goto  50 
add=add+gn(kk, ii-kk+1) *gd(kk) 

50  continue 

res=0 . 

do  55  jj=l,ll 

if(jj+ii-l  .gt.  ndof)goto  55 

res=res  +  gn(ii, j j) *gd( j j+ii-1) 

55  continue 

c 

c  add  to  existing  gf  which  already  contains  R 
c 

gf (ii) =gf0 (ii) * (pincr+tpincr) -res-add 
60  continue 
65  continue 
c 

c  impose  displacement  boundary  conditions 

c 

call  bndy ( ndo f,gsti00,gf, nbndry , ibndry , vbound ) 
c  if (icount.gt.l)call  bndy (ndof, gstif, gf, nbndry, ibndry, vpres) 

c 

c  solve  system  of  equations,  result  in  gf 

c 

call  solve (ndof , gstif , gf, 1, detm,detml) 
c 

c  through  line  69  copied  from  Tsai's  program 

c 

al=dss0 

a2=0.0 

a3=0.0 

do  147  ii=l,ndof 

a2=a2+ (gld(ii) +gf (ii) ) *gdis (ii) 

147  a3=a3+gf (ii) * (2 . 0*gld(ii) +gf (ii) ) 
dl2=a2*a2-al*a3 
c  write (6,*)  dl2,al,a2,a3 

c 
c 

if (dl2.lt. 0.0) then 
c 

c  deal  with  complex  roots  by  cutting  the  search 

c  radius  (dss)  in  half 

c 

do  2991  ii=l,ndof 
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2991  gd(ii) =gdOO (ii) 
iicut=iicut  +  1 
if (iicut .gt . 10) then 
write{6,3000) 
stop 
endif 
dss=dss/2 . 0 
goto  3 
endif 
iicut=0 

dpincl= ( -a2+dsqrt (dl2 ) ) /al 

dpinc2= (-a2-dsqrt (dl2) ) /al 

thetal=0 . 0 

theta2=0 . 0 

do  148  ii=l,ndof 

gldO (ii) =gld(ii) 

gld(ii) =gld(ii) +gf (ii) +dpincl*gdis (ii) 
gldl ( ii ) =gld0 ( ii ) +gf ( ii ) +dpinc2  *gdis ( ii ) 
thetal=thetal+gldO (ii) *gld(ii) 
theta2=theta2+gld0 (ii) *gldl (ii) 

148  continue 

c  write (6,*)  thetal , theta2 

thetl2=thetal*theta2 
if (thetl2 .gt.0.0)  go  to  149 
dpincr=dpincl 

if (theta2.gt.0.0)  call  chsign(gld,gldl,dpincr,dpinc2,ndof) 
go  to  150 

149  dpib=-a3/(a2*2.0) 
dpinl=dabs (dpib-dpincl) 
dpin2=dabs (dpib-dpinc2 ) 
dpincr=dpincl 

if (dpin2.lt. dpinl)  call  chsign(gld,gldl,dpincr,dpinc2,ndof ) 

150  pincr=pincr+dpincr 

69  continue 
c 

c  update  total  displacement  vector  gd 
c 

do  70  ii=l,ndof 

70  gd(ii) =gd(ii) +gld(ii) -gldO (ii) 
if (linear. eq.l) goto  80 

call  converge (ndof , neon, icount, tol, imax) 
c 

c  if  no  convergence  (ncon=0)  start  next  iteration 

c 

if (ncon.eq.O)goto  4 
80  continue 

if (neon. eq.l  .and.  ncount.le. nine) then 

call  postpr (icount, ncount, kcall, ndof ,pname, sname) 
i f ( ncount . eq . nine ) s  top 
ncount =ncount+l 
tpincr=tpincr+pincr 
goto  3 
endif 

3000  format (lx, 'More  than  10  consecutive  imaginary  roots') 
3010  format (/, lx, i2 ) 
return 
end 
c 
c 
c 
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subroutine  solve (ndof , band, rhs , ires , detm, detml ) 

C  . ; . 

c  solve  a  banded  symmetric  system  of  equations 

c  . . 

c 

implicit  double  precision  (a-h,o-2) 
dimension  band(2500, 11) ,rhs (2500) 
meqns=ndof-l 
if (ires . gt. 0) goto  90 
do  500  npiv=l,meqns 
c  print*, 'npiv=  ' ,npiv 

npivot=npiv+l 
lstsub=npiv+ll-l 
if (Istsub.gt .ndof )  lstsub=ndof 
do  400  nrow=npivot,  Istsub 
c  invert  rows  and  columns  for  row  factor 

ncol=nrow-npiv+l 

factor=band (npiv, ncol ) /band (npiv, 1) 
do  200  ncol=nrow, Istsub 
icol=ncol-nrow+l 
jcol=ncol-npiv+l 

200  band (nrow, icol ) =band (nrow, icol ) -f actor*band (npiv, j col ) 
400  rhs (nrow)  =rhs (nrow) -factor*rhs (npiv) 

500  continue 
detm=l . 0 
detml=0 . 0 
do  600  ii=l,ndof 

c  write  (*,*)  ' BAND ( ' , ii, ' 1) = ’ ,band(ii, 1) 

detml=detml+dlogl0 (dabs (band(ii,l) ) ) 

600  detm=detm*band (ii , 1) /dabs (band(ii, 1) ) 
go  to  101 

90  do  100  npiv=l,meqns 
npivot=npiv+l 
lstsub=npiv+ll-l 
if (Istsub. gt .ndof )  lstsub=ndof 
do  110  nrow=npivot,  Istsub . 
ncol=nrow-npiv+l 

factor=band (npiv, ncol) /band (npiv, 1) 

110  rhs (nrow) =rhs (nrow) -factor*rhs (npiv) 

100  continue 
c  back  substitution 

101  do  800  ijk=2,ndof 
np i v=ndo f - i j k+ 2 

rhs (npiv) =rhs (npiv) /band (npiv, 1) 

lstsub=npiv-ll+l 

if (Istsub. It. 1)  lstsub=l 

npivot=npiv-l 

do  700  jki=lstsub,npivot 

nrow=npivot- j ki+lstsub 

ncol=npiv-nrow+l 

f actor=band (nrow, ncol ) 

700  rhs (nrow) =rhs (nrow) -factor*rhs (npiv) 

800  continue 

rhs ( 1 ) =rhs ( 1 ) /band (1,1) 
return 
end 
c 
c 
c 

subroutine  chsign (gld, gldl , dpincr , dpinc2 , ndof ) 
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implicit  double  precision  (a-h,o-z) 
dimension  gld{2500) ,gldl (2500) 
do  100  i=l,ndof 
100  gld(i)=gldl (i) 
dpincr=dpinc2 
return 
end 
c 
c 
c 

subroutine  stiff (ielem, icount,ncount,kcall) 
c 

implicit  double  precision  {a-h,o-z) 
character*64  gname,fname 
common /chac/gname, fname 
c 

conffnon/ input /tol , table (250 ) , delem(250) ,vbound (2500 ) , distld, 

.  vconc (2500 ) , ey , enu, ht , el , e2 , gl2 , enul2 , enu21 , gl3 , g23 , pthick, 
rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 

.  nelem,nbndry,nbound(250, 6) , Idtyp.nconc, iconc (2500) , 

.  nplies,nforc,iforc(2500) ,nstres, istres (250) , ibndry (2500) , 
theta(20)  ,idload(250)  ,coord(251)  ,widt]i,nnod,pincr, eiter, ttpi 
c 

common/ elas /ae, de,  f e, he, e j , el, re, te, as , ds ,  f s 
c 

common/stf /stif (11,11) , elp (11) , eln.  (11, 11 ) , eld ( 11 ) 
c 

common /shp/dsf (7,11) 
c 

dimension  bmk (7,7) , bmnl (7,7) ,bmn2 (7,7) , 

.  gauss4(4) ,wt4(4) ,gauss7(7) ,wt7 (7) , g (7) , dsftr (11, 7) , 

pkt(7,7) ,pkn(7,7) ,pktd(7,ll) ,pknd(7,ll) ,gauss5(5) ,wt5(5) 

data  gauss4/0 . 86113 63 115d0 , 0 . 339981043 5d0 , -0 . 33 99 81043 5d0 , 
-0.8611363115d0/ 

data  wt4/0.3478548451d0,0.6521451548d0, 0. 6521451548d0, 
0.3478548451d0/ 

data  gauss5/0.9061798459d0,0.5384693101d0,0.0d0,-0.5384693101d0, 
-0.9061798459d0/ 

data  wt5/0.2369268851d0,0.4786286705d0,0.5688888889d0, 

0. 47862 86705d0, 0.23 692 68851d0/ 
data  gauss7/0.9491079123d0,0.7415311856d0, 0.4058451513d0, 

0 . OdO , -0 . 4058451513d0 , -0 . 7415311856d0 , -0 . 9491079123d0/ 
data  wt7/0.1294849662d0,0.2797053915d0,0.3818300505d0, 

0.4179591836d0,0.3818300505d0,0.2797053915d0,0.1294849662d0/ 

c 

c  initialize  stiffness  arrays  and  load  array 

c 

do  10  ii=l,ll 
elp (ii) =0 . 0 
do  10  jj=l,ll 
stif (ii, j j ) =0 . 0 
10  eln(ii, j j)=0.0 
c 

c  set  n\imber  of  gauss  points  for  interpolation 

c 

ngp=5 

if (ncount.eq.l  .and.  icount.eq. l)ngp=4 
if ( linear. eq.l)ngp=4 
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ekl=-4./(3.*ht**2) 

if  (isarch.eq.l)  pl=l./rad 

if (ncount .eq. 1  .and.  icount.eq.l 

.and. 

kcall.eq.l 

.  and 

isarch.eq.l) 

call  beaink(bmk, ekl,pl) 
if (ncount. eq.l  .and.  icount.eq.l 

.and. 

kcall.eq.l 

.  and 

isarch.eq. 0) 

call  sbeamk (bmk, ekl) 


c 

c  loop  over  gauss  points 

c 

do  100  ii=l,ngp 
if (ngp.eq.4)eta=gauss4 (ii) 
if (ngp.eq.5)eta=gauss5 (ii) 
if (ngp.eq.7)eta=gauss7 (ii) 
call  shape (eta, ielem, aa) 
c 

c  multiply  element  displacement  vector,  eld  (this  is  'q' 
c  in  the  thesis)  by  the  shape  function  matrix,  dsf,  to  get 
c  the  displacement  gradient  vector (d(s)  in  thesis,  q  here) 

c 

do  20  kk=l,7 
20  q(kk)=0.0 
c 

do  30  jj =1,7 
do  30  kk=l,ll 

30  q(j j)=q(j j)+dsf (j j,kk) *eld(kk) 


c 

c 

c  initialize  bmnl,  bmn2 
c 

do  35  kk=l,7 
do  35  jj=l,7 
bmnl ( j  j , kk) =0 . OdO 
35  bmn2 ( j j , kk) =0 . OdO 
c 

c  skip  bmnl  and  bmn2  comps  first  time  through 

c 

if(icount  ,eq.  1  .and.  ncovint  .eq.  l)goto  37 
c 

if (isarch.eq.l) call  beamnl (q, bmnl, ekl, pi) 
if (isarch.eq.l)call  beamn2 (q,bmn2, ekl,pl) 
if (isarch.eq.O)call  sbmnl (q, bmnl, ekl) 
if (isarch.eq.O)call  sbmn2 (q,bmn2,ekl) 

37  continue 


c 

c  transpose  the  shape  function  matrix 

c 

do  40  jj=l,7 
do  40  kk=l,ll 
40  dsftr (kk, j j ) =dsf ( j j ,kk) 
c 

c  create  element  independent  incremental  stiffness  array, 

c  pkt,  and  element  ind.  equilibrium  stiffness  array,  pkn 
c 

do  50  jj=l,7 
do  50  kk=l,7 

pkt ( j j ,kk) =bmk( j j ,kk) +bmnl ( j j ,kk) +bmn2 ( j j ,kk) 

50  pkn ( j j , kk) =bmk( j j , kk) +bmnl ( j j ,kk) /2 . +bmn2 ( j  j , kk) /3 . 
c 
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c  post-multiply  each  array  by  the  shape  function  matrix 

c 

do  60  jj=l,7 
do  60  kk=l,ll 
pktd ( j j , kk) =0 . 0 
pknd( j j , kk) =0 . 0 
do  60  11=1,7 

pktd(jj,kk)=pktd(jj,kk)  +  aa*pkt ( j j , 11) *dsf (11, kk) 

60  pknd( j j , kk) =pknd( j j ,kk)  +  aa*pkn( j j , 11) *dsf (11, kk) 
c 

c  Finally,  pre-multiply  these  new  arrays  by  the  transpose 

c  of  the  shape  function  matrix  to  get  the  element  incremental 

c  stiffness,  stif,  and  element  equilibrium  stiffness,  eln. 

c  Also  multiply  by  the  weighting  factor  for  this  particular 

c  gauss  point.  Note  that  these  arrays  are  zeroed  outside  the 
c  loop  over  the  gauss  points  since  they  accumulate  (integrate) 

c  data  over  all  the  gauss  points, 

c 

if (ngp.eq.4)wt=wt4 (ii) 
if (ngp.eq.5)wt=wt5 (ii) 
if (ngp.eq. 7) wt=wt7 (ii) 
do  70  jj=l,ll 
do  70  kk=l,ll 
do  70  11=1,7 

stif ( j j ,kk)=stif ( j j ,kk)+wt*width*dsftr (j j ,11) *pktd(ll,kk) 

eln ( j  j , kk) =eln ( j  j , kk) +wt *width*dsf tr ( j  j , 11 ) *pknd ( 11 , kk) 

70  continue 
100  continue 
c  write (6, 1010) ielem 

c  write (6, 1005) 

c  do  900  ii=l,ll 

c900  write(6,1000)  (stif (ii, j j ) , j j=l, 11) 

1000  format (9 (2x, dl2 . 5) ) 

1005  format (/,' stif ' ) 

1010  format (i4) 
return 
end 
c 
c 
c 

subroutine  shape (eta, ielem, aa) 
c 

implicit  double  precision  (a-h,o-z) 
c 
c 

common /shp/dSf (7,11) 
c 
c 

common/ input /tol , table (250 ) , delem(250) , vbound (2500 ) , distld, 
vconc (2500 ) , ey, enu, ht , el , e2 , gl2 , enul2 , enu21 , gl3 , g23 , pthick, 
rad, linear, isotro, isarch, ishape, inctyp,ninc, imax, 
nelem,nbndry,nbound(250, 6) , ldtyp,nconc, iconc (2500) , 

.  nplies,nforc,iforc(2500) ,nstres, istres (250) , ibndry (2500)  , 

.  theta(20) ,idload(250) ,coord(251) ,width,nnod,pincr, eiter,  ttpi 
c 

c  initialize  shape  function  matrix 

c 

do  10  ii=l,7 
do  10  jj=l,ll 
10  dsf (ii, j j ) =0 . 0 
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aa= (coord (ielem+1) -coord (ielem) ) *0 . 5 
c 

c  enter  values  into  dsf 

c  these  include  jacobian  terms 

c 
c 

c  Q1 , Q3 , Q2  and  derivatives 

c 

dsf (1,1) =0.5* (eta**2-eta) 
dsf (l,6)=1.0-eta**2 
dsf (1,7) =0.5* (eta**2+eta) 
dsf (2 , 1) = (eta-0 .5) /aa 
dsf (2,6)=-2.0*eta/aa 
dsf (2 , 7) = (eta+0 .5) /aa 
dsf (3,4)=0.25*(2.0-3.0*eta+eta**3) 
dsf (3 , 5) =0 . 25*aa* (1 . 0-eta-eta**2+eta**3 ) 
dsf (3, 10) =0.25* (2.0+3.0*eta-eta**3) 
dsf (3 , 11) =0 . 25*aa* (-1 . 0-eta+eta**2+eta**3 ) 
dsf (4, 4) =0.25* (-3.0+3.0*eta**2) /aa 
dsf (4, 5) =0.25* (-1.0-2.0*eta+3.0*eta**2) 
dsf (4,10)=0.25*(3.0-3.0*eta**2)/aa 
dsf (4, 11) =0.25* (-1.0+2 . 0*eta+3 .0*eta**2) 
dsf (5,4)=0.25*6.0*eta/aa**2 
dsf (5,5)=0.25*(-2.0+6.0*eta)/aa 
dsf (5,10)=-0.25*6.0*eta/aa**2 
dsf (5, 11) =0.25* (2.0+6.0*eta) /aa 
dsf (6, 2) =0.25* (2.0-3.0*eta+eta**3) 
dsf (6,3)=0.25*aa* (1.0-eta-eta**2+eta**3) 
dsf (6,8)=0.25*(2.0+3.0*eta-eta**3) 
dsf (6, 9)=0.25*aa* (-1 . 0-eta+eta**2+eta**3 ) 
dsf (7,2)=0.25*(-3.0+3.0*eta**2)/aa 
dsf (7, 3 ) =0 . 25* (-1 . 0-2 . 0*eta+3 . 0*eta**2 ) 
dsf (7, 8) =0.25* (3.0-3.0*eta**2) /aa 
dsf (7,9)=0.25*(-1.0+2.0*eta+3.0*eta**2) 
c 

c  temporary  print 

c 

c  do  100  ii=l,7 

c  100  write(6, 1000)  (dsf (ii, j j ) , j j=l, 11) 
c  1000  format  (9 (2x,dl2 .5) ) 

return 
end 
c 
c 
c 

subroutine  postpr (icount , ncovmt , kcall, ndof , pname, sname) 
c 

implicit  double  precision  (a-h,o-z) 
c 

character*64  gname,fname 
c 
c 

common /chac/gname, fname 
c 

common/elas/ae , de, f e , he, e j , el , re, te, as , ds , f s 
c 

common /input/ to 1 , table (250 ) , delem(250) , vbound (2500 ) , distld, 

.  vconc  (2500)  ,  ey,  enu, ht ,  el ,  e2 ,  gl2 ,  enul2 ,  enu21 ,  gl3 ,  g23 , pthiclc, 
rad, linear, isotro, isarch, ishape, inctyp,ninc, imax. 
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nelem,nbndry,nbound(250,6) , ldtyp,nconc, iconc (2500) , 
nplies,nforc,iforc(2500) ,nstres, istres (250) , ibndry (2500)  , 
theta(20) , idload(250) ,coord(251) , width, nnod,pincr, eiter, ttpi 

common/ stf/stif (11,11) , elp (11) , eln (11 , 11) , eld (11) 


c 

0) 


c 

c 

c 

c 


common/proc/gstif (2500,11) ,gn (2500, 11) , gf (2500) ,gd(2500) ,vperm(250 
vpres(2500) 

dimension  vforc(2500) ,xcoord(251) ,ycoord(251) 
pi=3. 14159 

if  ishape=l  print  global  x,y  coords  to  file  bshape 
if  ishape=2  figure  out  symmetric  coords  as  well 


if  (isarch.eq. 0)  then 
go  to  423 
end  if 

i f ( ishape . eq . 1 . or . i shape . eq . 2 . and . ncount . eq . 1 . and . isarch . eq . 1 ) then 
open (8, file=sname,status='new' ) 
do  1  ii=l,nnod 
if ( i shape. eq. 2) then 

xcoord(nnod-l+ii) =rad*cos (pi/2 . 0-coord (ii) /rad) 
xcoord (nnod+l-ii ) =-rad*cos (pi/2 . 0-coord ( ii ) /rad) 
ycoord(nnod-l+ii) =rad*sin(pi/2 . 0-coord (ii) /rad) 
ycoord (nnod+l-ii ) =ycoord (nnod-l+ii ) 
else 

xcoord(ii)=-rad*cos(pi/2.0+coord(ii) /rad- 
coord  (nnod)  /  (2*rad)  ) 

ycoord (ii)=rad*sin (pi/2. 0+coord( 11) /rad-coord (nnod) / (2*rad) ) 
endif 

1  continue 

do  100  ii=l,nnod 

100  write (8, 2000)xcoord(ii) ,ycoord(ii) 

if (ishape . eq. 2 ) then 
do  110  ii=2,nnod 

110  write (8,2000) xcoord (nnod+ii-1 ) , ycoord (nnod+ii-1 ) 

endif 

write (8 , 2010 ) 


c 

c 

c 


endif 

global  displacements  for  straight  beams 

423  if (ishape. eq.l .and. ncount .eq.l. and. isarch.eq. 0) then 
open (8, file=sname, status='new' ) 
do  2  ii=l,nnod 

xcoord (ii)=  coord (nnod) -coord (ii) 
ycoord (ii) =0 . 0 

2  write ( 8, 2 000) xcoord (ii) , ycoord (ii) 

write (8,2010) 

endif 


c 

c 

c 


print  out  global  displacements 

write (6, 1000) 
write  ( 6 , 1010 ) 

write ( 6 , 1020) ncount , icount 
write (6, 1030) 
do  90  ii=0,nnod-l 
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x,y  for  arches 


if ( ishape . eq . 1 . and . isarch . eq . 1 ) then 

i f (ishape . eq . 1 . and . isarch . eq . 1 . and .ii.eq.0)write(8,2020) ncount 
xcoord(ii+l) =- (rad-gd (6*ii+4) ) * 

cos (pi /2 . 0+coord ( ii+1 ) /rad-coord (nnod) / ( 2  *rad)  + 
gd(6*ii+l) /rad) 

ycoord(ii+l) = (rad-gd(6*ii+4) ) * 

sin (pi/2 .0+coord(ii+l) /rad-coord (nnod) / (2*rad)+ 
gd(6*ii+l) /rad) 

write (8, 2000)xcoord(ii+l) .ycoord (ii+1) 

endif 

i f ( i shape . eq . 2 . and . isarch . eq . 1 ) then 
if (ii.eq. 0) write (8, 2020) ncount 

xcoord(nnod+ii) = (rad-gd(6*ii+4) ) * 

cos (pi/2. 0-coord(ii+l) /rad  +gd(6*ii+l) /rad) 
xcoord (nnod-ii ) =-xcoord (nnod+ii ) 
ycoord(nnod+ii) = (rad-gd(6*ii+4) ) * 

sin(pi/2 . 0-coord (ii+1) /rad+gd(6*ii+l) /rad) 
ycoord (nnod-ii ) =ycoord (nnod+ii) 

endif 

x,y  for  straight  beams 

if ( ishape. eq. 1 .and. isarch. eq. 0) then 

if (ishape . eq. 1 .and. isarch. eq. 0 .and. ii .eq. 0)  write ( 8 , 2 020) ncount 
xcoord(ii+l) =coord(nnod) -coord(ii+l) -gd(6*ii+l) 
ycoord(ii+l) =-gd(6*ii+4) 
write (8, 2000) xcoord (ii+1) , ycoord (ii+1) 

endif 

write(6, 1040) ii+1, (gd(6*ii+jj) , jj=l,5) 
write (6,1050)gd(6*ii+6) 
if ( ishape . eq. 2 . and. isarch. eq. 1) then 
do  95  ii=l, 2*nnod-l 
write (8, 2000) xcoord (ii) , ycoord (ii) 

endif 

if (ishape .ge . 1) write (8,2010) 

compute  equivalent  forces  requested 

do  5  ii=l,ndof 
gf (ii) =0 . OdO 
do  5  jj=l,ll 
gstif (ii , j j ) =0 . OdO 
gn(ii, jj)=0.0d0 

loop  over  all  elements  for  stiffness  and  forces 

do  30  ielem=l,nelem 

do  10  ii=l,ll 

eld(ii) =gd(ii+ (ielem-1) *6) 

call  stiff (ieiem, icount, ncount, kcall) 

Assemble  global  stiffness  array,  gstif,  global  equilibrium 
stiffness,  gn,  in  banded  form.  Half-bandwidth=ll .  Also 
assemble  global  force  vector,  gf. 


nr= (ielem-1) *6  +  1 


do  30  jj=0,10 

gf (nr+j j ) =gf (nr+j j ) +elp( j j+1) 
do  30  kk=l,ll-jj 

gstif (nr+j  j , kk) =gstif (nr+j  j , kk) +stif ( j  j  +1 , kk+ j  j ) 
if  (linear. eq.  Dgoto  30 

gn (nr+j j ,kk) =gn(nr+j j , kk) +eln( j j+l,kk+j j ) 

30  continue 
c 

c  calculate  the  residual  force  vector  for  nonlinear 

c  analysis.  - [gn] * {gd}+R=- [k+nl/2+n2/3 ] * {q}+R=gf 
c 

do  60  jj=l,nforc 
ii=iforc ( j j ) 
add=0 . 

do  50  kk=l,ii-l 
if(ii-kk+l  .gt.  11) goto  50 
add=add+gn (kk, ii-kk+1) *gd(kk) 

50  continue 

res=0 . 

do  55  11=1,11 

if(ll+ii-l  .gt.  ndof)goto  55 
res=res  +  gn(ii, 11) *gd(ll+ii-l) 

55  continue 

c 

c  compute  nodal  force 
c 

vf ore ( j  j ) =res+add 
60  continue 
c 

c  print  nodal  forces  and  create  plot  file 
c 

open (7, file=pname, status= 'new' ) 
if (ncount . eq. 1) write (7,1065)0.0,0.0,0.0 
write (6, 1060) 
do  70  ii=l,nforc 

write (7, 1065) gd(iforc (ii) ) ,gd(iforc (ii) -3 ) , vforc (ii) 

70  write (6, 1070) iforc(ii) , vforc (ii) 

1000  format (/) 

1010  format (lx, 'Results  of  nonlinear  analysis') 

1020  format (lx, ' increment= ', i3 , '  '  iteration= ' , i3 ) 

1030  format (lx,  'Node' ,7x, 'V , 12x,  'Psi-s' , 9x, 'Psi-ss' ,9x, 'W ,9x, 'W-s 
1040  format (lx, i4 , 5 (2x, dl2 . 5) ) 

1050  format (lx, 'Midnode  v; ' ,3x,dl2 .5) 

1060  format (lx, /,' Equivalent  nodal  forces:') 

1065  format (lx, dl2 . 5, 2x, dl2 . 5 , 2x,dl2 .5) 

1070  format ( lx, ' DOF  no : ' , i4 , 2x, ' Force : ' , lx, dl2 . 5 ) 

2000  format (lx, fl2 . 5 , 2x, fl2 . 5) 

2010  format (//) 

2020  format (/, lx, i4) 
return 
end 
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